{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 135,
   "id": "6a4298c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from matplotlib.lines import Line2D\n",
    "from linearmodels.panel import PanelOLS\n",
    "import statsmodels.api as sm\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5cfa964f",
   "metadata": {},
   "source": [
    "## Create simulation data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 136,
   "id": "cd338178",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set seed\n",
    "np.random.seed(20200403)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 137,
   "id": "a11c9789",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1,000 firms (25 per state), 40 states, 4 groups (250 per groups), 30 years\n",
    "\n",
    "df = pd.DataFrame({\"state\" : np.repeat(np.arange(1, 41), 25),\n",
    "                   \"firms\" : np.random.uniform(0, 5, size = 1000)})\n",
    "\n",
    "df = pd.concat([df] * 30, ignore_index = True).sort_values(['state', 'firms'])\n",
    "\n",
    "df['year'] = list(range(1980, 2010)) * 1000\n",
    "\n",
    "df['n'] = list(range(1, 31)) * 1000\n",
    "\n",
    "df['id'] = np.repeat(np.arange(1, 1001), 30)\n",
    "\n",
    "df.reset_index(drop = True, inplace = True)\n",
    "\n",
    "df['group'] = np.repeat(np.arange(1, 5), len(df)/4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 139,
   "id": "0bade3fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add 250 firms treated every period with the treatment effect still 7 on average\n",
    "# Cohort years 1986, 1992, 1998, 2004\n",
    "\n",
    "def create_treat_date(row):\n",
    "    if row['group'] == 1:\n",
    "        return 1986\n",
    "    elif row['group'] == 2:\n",
    "        return 1992\n",
    "    elif row['group'] == 3:\n",
    "        return 1998\n",
    "    else:\n",
    "        return 2004\n",
    "    \n",
    "df['treat_date'] = df.apply(create_treat_date, axis = 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 140,
   "id": "823aa3ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_treat(row):\n",
    "    if (row['group'] == 1 and row['year'] >= 1986 or  \n",
    "        row['group'] == 2 and row['year'] >= 1992 or \n",
    "        row['group'] == 3 and row['year'] >= 1998 or \n",
    "        row['group'] == 4 and row['year'] >= 2004):\n",
    "        return 1\n",
    "    else:\n",
    "        return 0\n",
    "    \n",
    "df['treat'] = df.apply(create_treat, axis = 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 142,
   "id": "4850e1cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Data generating process \n",
    "\n",
    "def create_te(row):\n",
    "    if row['group'] == 1:\n",
    "        return np.random.normal(10, 0.2**2)\n",
    "    elif row['group'] == 2:\n",
    "        return np.random.normal(8, 0.2**2)\n",
    "    elif row['group'] == 3:\n",
    "        return np.random.normal(6, 0.2**2)\n",
    "    else:\n",
    "        return np.random.normal(4, 0.2**2)  \n",
    "    \n",
    "df['te'] = df.apply(create_te, axis = 1)\n",
    "df['e'] = np.random.normal(0, 0.5**2, size = 30000)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b841f557",
   "metadata": {},
   "source": [
    "* DGP: heterogeneous versus constant (but always across group heterogeneity)\n",
    "* Cumulative treatment effect is te x (year - t_g + 1) -- Dynamic treatment effects over time for each group.\n",
    "* How does (year - treat_date + 1) create dynamic ATT?  Assume treat_date is 1992 and it is year 2000. Then, te=8 x (2000 - 1992 + 1) = 8 x (9) = 72. Group 2's TE rises from an 8 up to 72 in the t+8 year."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 221,
   "id": "f4334074",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Constant treatment effects.  Notice, the treatment effect is constant. \n",
    "df['y2'] = (df['firms'] + df['n'] + df['te']*df['treat'] + df['e'])\n",
    "# Data generating process with heterogeneity over time\n",
    "df['y'] = (df['firms'] + df['n'] + df['treat']*df['te']*(\n",
    "    df['year'] - df['treat_date'] + 1) + df['e'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 144,
   "id": "f88167d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "df['individual'] = df['id']\n",
    "df['time'] = df['year']\n",
    "df = df.set_index(['individual', 'time'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7086eeb",
   "metadata": {},
   "source": [
    "* For group 1, the ATT in 1986 is 10\n",
    "* For group 1, the ATT in 1987 is 20\n",
    "* For group 1, the ATT in 1988 is 30 and so on\n",
    "* This is what we mean by \"dynamic treatment effects\" or \"heterogeneity over time\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "942a4219",
   "metadata": {},
   "source": [
    "## Estimation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "id": "6ee91f7c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                          PanelOLS Estimation Summary                           \n",
      "================================================================================\n",
      "Dep. Variable:                     y2   R-squared:                        0.9943\n",
      "Estimator:                   PanelOLS   R-squared (Between):              0.6274\n",
      "No. Observations:               30000   R-squared (Within):               0.9943\n",
      "Date:                Tue, Feb 15 2022   R-squared (Overall):              0.9688\n",
      "Time:                        21:45:57   Log-likelihood                -3.782e+04\n",
      "Cov. Estimator:                Robust                                           \n",
      "                                        F-statistic:                   1.682e+05\n",
      "Entities:                        1000   P-value                           0.0000\n",
      "Avg Obs:                       30.000   Distribution:                F(30,28970)\n",
      "Min Obs:                       30.000                                           \n",
      "Max Obs:                       30.000   F-statistic (robust):          1.581e+05\n",
      "                                        P-value                           0.0000\n",
      "Time periods:                      30   Distribution:                F(30,28970)\n",
      "Avg Obs:                      1000.00                                           \n",
      "Min Obs:                      1000.00                                           \n",
      "Max Obs:                      1000.00                                           \n",
      "                                                                                \n",
      "                                Parameter Estimates                                \n",
      "===================================================================================\n",
      "                 Parameter  Std. Err.     T-stat    P-value    Lower CI    Upper CI\n",
      "-----------------------------------------------------------------------------------\n",
      "Intercept           3.4148     0.0391     87.364     0.0000      3.3382      3.4915\n",
      "C(year)[T.1981]     1.0052     0.0552     18.208     0.0000      0.8970      1.1134\n",
      "C(year)[T.1982]     2.0119     0.0552     36.442     0.0000      1.9037      2.1201\n",
      "C(year)[T.1983]     2.9940     0.0554     54.063     0.0000      2.8855      3.1026\n",
      "C(year)[T.1984]     3.9869     0.0554     71.998     0.0000      3.8783      4.0954\n",
      "C(year)[T.1985]     4.9881     0.0553     90.271     0.0000      4.8798      5.0964\n",
      "C(year)[T.1986]     6.7422     0.0430     156.95     0.0000      6.6580      6.8264\n",
      "C(year)[T.1987]     7.7450     0.0429     180.49     0.0000      7.6609      7.8291\n",
      "C(year)[T.1988]     8.7413     0.0430     203.11     0.0000      8.6569      8.8256\n",
      "C(year)[T.1989]     9.7367     0.0429     227.10     0.0000      9.6526      9.8207\n",
      "C(year)[T.1990]     10.744     0.0430     249.75     0.0000      10.660      10.828\n",
      "C(year)[T.1991]     11.754     0.0429     273.69     0.0000      11.670      11.838\n",
      "C(year)[T.1992]     13.001     0.0403     322.50     0.0000      12.922      13.080\n",
      "C(year)[T.1993]     14.000     0.0403     347.72     0.0000      13.921      14.078\n",
      "C(year)[T.1994]     15.009     0.0403     372.33     0.0000      14.930      15.088\n",
      "C(year)[T.1995]     16.000     0.0403     397.37     0.0000      15.921      16.078\n",
      "C(year)[T.1996]     17.006     0.0402     422.81     0.0000      16.927      17.084\n",
      "C(year)[T.1997]     18.008     0.0402     447.88     0.0000      17.929      18.087\n",
      "C(year)[T.1998]     18.743     0.0438     427.85     0.0000      18.657      18.829\n",
      "C(year)[T.1999]     19.754     0.0437     452.11     0.0000      19.668      19.839\n",
      "C(year)[T.2000]     20.749     0.0437     474.59     0.0000      20.663      20.835\n",
      "C(year)[T.2001]     21.751     0.0436     498.69     0.0000      21.665      21.836\n",
      "C(year)[T.2002]     22.741     0.0438     519.42     0.0000      22.655      22.827\n",
      "C(year)[T.2003]     23.756     0.0438     541.84     0.0000      23.670      23.841\n",
      "C(year)[T.2004]     23.999     0.0559     429.47     0.0000      23.890      24.109\n",
      "C(year)[T.2005]     25.000     0.0560     446.48     0.0000      24.890      25.110\n",
      "C(year)[T.2006]     26.011     0.0559     465.49     0.0000      25.902      26.121\n",
      "C(year)[T.2007]     26.992     0.0561     480.71     0.0000      26.882      27.102\n",
      "C(year)[T.2008]     28.014     0.0561     499.55     0.0000      27.905      28.124\n",
      "C(year)[T.2009]     29.002     0.0559     519.10     0.0000      28.893      29.112\n",
      "treat               6.9948     0.0193     362.46     0.0000      6.9570      7.0327\n",
      "===================================================================================\n",
      "\n",
      "F-test for Poolability: 120.98\n",
      "P-value: 0.0000\n",
      "Distribution: F(999,28970)\n",
      "\n",
      "Included effects: Entity\n"
     ]
    }
   ],
   "source": [
    "# Estimation using TWFE - constant treatment effects\n",
    "constant_te = PanelOLS.from_formula('y2 ~ 1 + treat + C(year) + EntityEffects', \n",
    "                                    data = df).fit(cov_type = 'robust')\n",
    "print(constant_te)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 146,
   "id": "d5cdf451",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                          PanelOLS Estimation Summary                           \n",
      "================================================================================\n",
      "Dep. Variable:                      y   R-squared:                        0.7229\n",
      "Estimator:                   PanelOLS   R-squared (Between):             -0.0784\n",
      "No. Observations:               30000   R-squared (Within):               0.7229\n",
      "Date:                Tue, Feb 15 2022   R-squared (Overall):              0.4704\n",
      "Time:                        21:45:58   Log-likelihood                -1.437e+05\n",
      "Cov. Estimator:                Robust                                           \n",
      "                                        F-statistic:                      2519.8\n",
      "Entities:                        1000   P-value                           0.0000\n",
      "Avg Obs:                       30.000   Distribution:                F(30,28970)\n",
      "Min Obs:                       30.000                                           \n",
      "Max Obs:                       30.000   F-statistic (robust):             2084.7\n",
      "                                        P-value                           0.0000\n",
      "Time periods:                      30   Distribution:                F(30,28970)\n",
      "Avg Obs:                      1000.00                                           \n",
      "Min Obs:                      1000.00                                           \n",
      "Max Obs:                      1000.00                                           \n",
      "                                                                                \n",
      "                                Parameter Estimates                                \n",
      "===================================================================================\n",
      "                 Parameter  Std. Err.     T-stat    P-value    Lower CI    Upper CI\n",
      "-----------------------------------------------------------------------------------\n",
      "Intercept           3.4148     1.2512     2.7292     0.0064      0.9624      5.8673\n",
      "C(year)[T.1981]     1.0052     1.7695     0.5681     0.5700     -2.4630      4.4734\n",
      "C(year)[T.1982]     2.0119     1.7695     1.1370     0.2555     -1.4564      5.4803\n",
      "C(year)[T.1983]     2.9940     1.7696     1.6919     0.0907     -0.4745      6.4626\n",
      "C(year)[T.1984]     3.9869     1.7696     2.2530     0.0243      0.5184      7.4553\n",
      "C(year)[T.1985]     4.9881     1.7695     2.8189     0.0048      1.5198      8.4564\n",
      "C(year)[T.1986]     10.167     1.6105     6.3131     0.0000      7.0104      13.324\n",
      "C(year)[T.1987]     13.670     1.5430     8.8590     0.0000      10.645      16.694\n",
      "C(year)[T.1988]     17.166     1.4843     11.565     0.0000      14.257      20.075\n",
      "C(year)[T.1989]     20.661     1.4360     14.388     0.0000      17.847      23.476\n",
      "C(year)[T.1990]     24.170     1.3990     17.277     0.0000      21.428      26.912\n",
      "C(year)[T.1991]     27.675     1.3743     20.138     0.0000      24.982      30.369\n",
      "C(year)[T.1992]     34.840     1.3230     26.334     0.0000      32.247      37.433\n",
      "C(year)[T.1993]     40.352     1.2962     31.131     0.0000      37.812      42.893\n",
      "C(year)[T.1994]     45.861     1.2859     35.666     0.0000      43.340      48.381\n",
      "C(year)[T.1995]     51.341     1.2923     39.728     0.0000      48.808      53.873\n",
      "C(year)[T.1996]     56.867     1.3155     43.227     0.0000      54.288      59.445\n",
      "C(year)[T.1997]     62.358     1.3539     46.056     0.0000      59.704      65.011\n",
      "C(year)[T.1998]     71.018     1.3982     50.793     0.0000      68.277      73.758\n",
      "C(year)[T.1999]     78.039     1.4413     54.145     0.0000      75.214      80.864\n",
      "C(year)[T.2000]     85.018     1.4942     56.899     0.0000      82.089      87.947\n",
      "C(year)[T.2001]     92.025     1.5567     59.114     0.0000      88.973      95.076\n",
      "C(year)[T.2002]     99.031     1.6281     60.827     0.0000      95.840      102.22\n",
      "C(year)[T.2003]     106.05     1.7064     62.146     0.0000      102.70      109.39\n",
      "C(year)[T.2004]     115.69     1.7176     67.353     0.0000      112.32      119.05\n",
      "C(year)[T.2005]     123.72     1.7622     70.207     0.0000      120.26      127.17\n",
      "C(year)[T.2006]     131.68     1.8075     72.848     0.0000      128.13      135.22\n",
      "C(year)[T.2007]     139.70     1.8559     75.272     0.0000      136.06      143.33\n",
      "C(year)[T.2008]     147.75     1.9057     77.529     0.0000      144.01      151.48\n",
      "C(year)[T.2009]     155.70     1.9549     79.645     0.0000      151.86      159.53\n",
      "treat              -6.7043     0.6662    -10.063     0.0000     -8.0100     -5.3985\n",
      "===================================================================================\n",
      "\n",
      "F-test for Poolability: 32.683\n",
      "P-value: 0.0000\n",
      "Distribution: F(999,28970)\n",
      "\n",
      "Included effects: Entity\n"
     ]
    }
   ],
   "source": [
    "# Estimation using TWFE - heterogenous treatment effects over time\n",
    "het_te = PanelOLS.from_formula('y ~ 1 + treat + C(year) + EntityEffects', \n",
    "                               data = df).fit(cov_type = 'robust')\n",
    "print(het_te)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35e23b3d",
   "metadata": {},
   "source": [
    "### Sun and Abraham event study commentary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 147,
   "id": "e13fcbbe",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Leads and lags\n",
    "df['time_til'] = df['year'] - df['treat_date']\n",
    "df['cons'] = 1\n",
    "df = pd.concat([df, pd.get_dummies(df['time_til'], prefix = \"dd\")], axis = 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 148,
   "id": "31bbf281",
   "metadata": {},
   "outputs": [],
   "source": [
    "dd = [dd for dd in df.columns if dd.startswith(\"dd_\")]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ba1f702",
   "metadata": {},
   "source": [
    "##### Event study with heterogeneity, dropping two leads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 150,
   "id": "82dc6aa8",
   "metadata": {},
   "outputs": [],
   "source": [
    "leads = {'dd_-4', 'dd_-1'}\n",
    "ind_cols = [ind for ind in dd if ind not in leads]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 167,
   "id": "73a0fb5f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6cAAAJXCAYAAABv+SLdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABAHElEQVR4nO3dfZRsa10f+O/TXCWmulQIL+IBvMQuY8BRZI7ojJ0RYxTSaxRNMIXOGIzHkJmoK3icWcCYxLcQCEk0xiXOIlcNmhgKRYWwWgVRYtoY4VwEeU9XRODaBK5RY3VNRKGe+aPrNE3fOuee09W7du/uz2etXl29d9WvfqdO7ar61rP3s0utNQAAANCmtbYbAAAAAOEUAACA1gmnAAAAtE44BQAAoHXCKQAAAK0TTgEAAGjdHW03cNRDHvKQeuedd7bdBgAAAA24++67f7fW+tBF685UOL3zzjtz7dq1ttsAAACgAaWU99xond16AQAAaJ1wCgAAQOuEUwAAAFonnAIAANA64RQAAIDWCacAAAC0TjgFAACgdcIpAAAArRNOAQAAaJ1wCgAAQOuEUwAAAFonnAIAANA64RQAAIDWCacAAAC0TjgFAACgdcIpAAAArRNOAQAAaJ1wCgAAQOuEUwAAAFonnAIAANC6O9pugPNvMplkNBpld3c3g8Egw+Ew/X6/7bYAAIAzRDilUTs7O9na2spsNst0Ok2v18vVq1ezvb2dzc3NttsDAADOCLv10pjJZJKtra1MJpNMp9MkyXQ6PVy+v7/fcocAAMBZIZzSmNFolNlstnDdbDbLaDRacUcAAMBZJZzSmN3d3cMR0+Om02nG4/GKOwIAAM4q4ZTGDAaD9Hq9het6vV42NjZW3BEAAHBWCac0ZjgcZm1t8VNsbW0tw+FwxR0BAABnlXBKY/r9fra3t9Pv9w9HUHu93uHy9fX1ljsEAADOCqeSoVGbm5vZ29vLaDTKeDzOxsZGhsOhYAoAAHwM4ZTGra+v58qVK223AQAAnGF26wUAAKB1wikAAACtE04BAABonXAKAABA64RTAAAAWiecAgAA0DrhFAAAgNYJpwAAALROOAUAAKB1wikAAACtE04BAABonXAKAABA64RTAAAAWiecAgAA0DrhFAAAgNYJpwAAALROOAUAAKB1txxOSyl/qpTy+lLKm0spbyulfNd8+YNLKa8ppezOfz/oyG2eW0oZl1LeVUp5chP/AAAAALrvdkZOP5TkL9ZaPyfJ45M8pZTyBUmek+S1tdZBktfO/04p5bFJnp7kcUmekuRFpZQHnGLvAAAAnBO3HE7rgf35nx83/6lJnprkJfPlL0nylfPLT03y0lrrh2qt704yTvLE02gaAACA8+W2jjktpTyglPKmJB9M8ppa668neXit9f1JMv/9sPnVLyV535Gb3zNfdrzmM0sp10op1+69994T/BMAAADoutsKp7XWj9RaH5/kkUmeWEr5rJtcvSwqsaDmi2utl2utlx/60IfeTjsAAACcEyearbfW+gdJXpeDY0k/UEp5RJLMf39wfrV7kjzqyM0emWTvpI0CAABwft3ObL0PLaV88vzyJyT5S0nemeSVSZ4xv9ozkrxifvmVSZ5eSnlgKeUxSQZJXn9KfQMAAHCO3HEb131EkpfMZ9xdS/KyWuurSim/luRlpZQrSd6b5KuTpNb6tlLKy5K8PcmHk3xTrfUjp9s+AAAA50Gp9T6Hgbbm8uXL9dq1a223AQAAQANKKXfXWi8vWneiY04BAADgNAmnAAAAtE44BQAAoHXCKQAAAK0TTgEAAGidcAoAAEDrhFMAAABaJ5wCAADQOuEUAACA1gmnAAAAtE44BQAAoHXCKQAAAK0TTgEAAGidcAoAAEDrhFMAAABaJ5wCAADQOuEUAACA1gmnAAAAtE44BQAAoHXCKQAAAK0TTgEAAGidcAoAAEDrhFMAAABaJ5wCAADQOuEUAACA1gmnAAAAtE44BQAAoHXCKQAAAK0TTgEAAGidcAoAAEDrhFMAAABaJ5wCAADQOuEUAACA1gmnAAAAtE44BQAAoHXCKQAAAK0TTgEAAGidcAoAAEDrhFMAAABaJ5wCAADQOuEUAACA1gmnAAAAtE44BQAAoHXCKQAAAK0TTgEAAGidcAoAAEDrhFMAAABaJ5wCAADQOuEUAACA1gmnAAAAtE44BQAAoHXCKQAAAK27o+0GbtdkMsloNMru7m4Gg0GGw2H6/X7bbQG3wPYLAMCNlFpr2z0cunz5cr127doN1+/s7GRrayuz2SzT6TS9Xi9ra2vZ3t7O5ubmCjsFbpftFwCAUsrdtdbLC9d1JZxOJpNcunQpk8nkPuv6/X729vayvr7edIvACdh+AQBIbh5OO3PM6Wg0ymw2W7huNptlNBqtuCPgVtl+AQC4P7ccTkspjyql/HIp5R2llLeVUv7OfPl3llJ+p5TypvnP1pHbPLeUMi6lvKuU8uRlGt3d3c10Ol24bjqdZjweL1MeaJDtFwCA+3M7EyJ9OMm31VrfWErpJ7m7lPKa+brvq7X+k6NXLqU8NsnTkzwuyacm+cVSymfUWj9ykkYHg0F6vd7CD7i9Xi8bGxsnKQusgO0XAID7c8sjp7XW99da3zi/PEnyjiSXbnKTpyZ5aa31Q7XWdycZJ3niSRsdDodZW1vc7traWobD4UlLAw2z/QIAcH9OdMxpKeXOJJ+b5Nfni765lPKbpZQfKaU8aL7sUpL3HbnZPbl5mL2pfr+f7e3t9Pv99Hq9JAcjLteXm0wFzi7bLwAA9+e2z3NaSllP8vIkz6q1/mEp5YeSfE+SOv/9T5N8Q5Ky4Ob3mRq4lPLMJM9Mkkc/+tE3ve/Nzc3s7e1lNBplPB5nY2Mjw+HQB1voANsvAAA3c1unkimlfFySVyX5hVrr9y5Yf2eSV9VaP6uU8twkqbU+f77uF5J8Z631125U//7OcwoAAEB3ncqpZEopJckPJ3nH0WBaSnnEkat9VZK3zi+/MsnTSykPLKU8Jskgyetvt3kAAADOv9vZrfcLk3xdkreUUt40X/b/JPmaUsrjc7DL7m8n+VtJUmt9WynlZUnenoOZfr/ppDP1AgAAcL7dcjitte5k8XGk2ze5zfOSPO8EfQEAAHCBnGi2XgAAADhNwikAAACtE04BAABonXAKAABA64RTAAAAWiecAgAA0DrhFAAAgNYJpwAAALROOAUAAKB1wikAAACtE04BAABonXAKAABA64RTAAAAWiecAgAA0DrhFAAAgNYJpwAAALROOAUAAKB1wikAAACtE04BAABonXAKAABA64RTAAAAWiecAgAA0DrhFAAAgNYJpwAAALTujrYbOO8mk0lGo1F2d3czGAwyHA7T7/fbbgsAAOBMEU4btLOzk62trcxms0yn0/R6vVy9ejXb29vZ3Nxsuz0AAIAzw269DZlMJtna2spkMsl0Ok2STKfTw+X7+/stdwgAAHB2CKcNGY1Gmc1mC9fNZrOMRqMVdwQAAHB2CacN2d3dPRwxPW46nWY8Hq+4IwAAgLNLOG3IYDBIr9dbuK7X62VjY2PFHQEAAJxdwmlDhsNh1tYWP7xra2sZDocr7ggAAODsEk4b0u/3s729nX6/fziC2uv1Dpevr6+33CEAAMDZ4VQyDdrc3Mze3l5Go1HG43E2NjYyHA4FUwAAgGOE04atr6/nypUrbbcBAABwptmtFwAAgNYJpwAAALROOAUAAKB1wikAAACtE04BAABonXAKAABA64RTAAAAWiecAgAA0DrhFAAAgNbd0XYDnMxkMsloNMru7m4Gg0GGw2H6/X7bbQEAAJyIcNpBOzs72draymw2y3Q6Ta/Xy9WrV7O9vZ3Nzc222wMAALhtduvtmMlkkq2trUwmk0yn0yTJdDo9XL6/v99yhwAAALdPOO2Y0WiU2Wy2cN1sNstoNFpxRwAAAMsTTjtmd3f3cMT0uOl0mvF4vOKOAAAAliecdsxgMEiv11u4rtfrZWNjY8UdAQAALE847ZjhcJi1tcX/bWtraxkOhyvuCAAAYHnCacf0+/1sb2+n3+8fjqD2er3D5evr6y13CAAAcPucSqaDNjc3s7e3l9FolPF4nI2NjQyHQ8EUAADoLOG0o9bX13PlypW22wAAADgVdusFAACgdcIpAAAArRNOAQAAaJ1wCgAAQOtuOZyWUh5VSvnlUso7SilvK6X8nfnyB5dSXlNK2Z3/ftCR2zy3lDIupbyrlPLkJv4BAAAAdN/tjJx+OMm31Vr/fJIvSPJNpZTHJnlOktfWWgdJXjv/O/N1T0/yuCRPSfKiUsoDTrN5AAAAzodbPpVMrfX9Sd4/vzwppbwjyaUkT03ypPnVXpLkdUmePV/+0lrrh5K8u5QyTvLEJL92Ws2fpslkktFolN3d3QwGgwyHw/T7/bbbAgAAuBBOdJ7TUsqdST43ya8nefg8uKbW+v5SysPmV7uU5D8eudk982Vnzs7OTra2tjKbzTKdTtPr9XL16tVsb29nc3Oz7fYAAADOvdsOp6WU9SQvT/KsWusfllJueNUFy+qCes9M8swkefSjH3277SxtMplka2srk8nkcNl0Ok2SbG1tZW9vL+vr6yvvqw1GjwEAgLbc1my9pZSPy0Ew/de11p+eL/5AKeUR8/WPSPLB+fJ7kjzqyM0fmWTveM1a64trrZdrrZcf+tCH3m7/SxuNRpnNZgvXzWazjEajFXfUjp2dnVy6dCnPetaz8sIXvjDPetazcunSpezs7LTdGgAAcAHczmy9JckPJ3lHrfV7j6x6ZZJnzC8/I8krjix/einlgaWUxyQZJHn98i2frt3d3cOR0uOm02nG4/GKO1q9o6PH1x+L6XR6uHx/f7/lDgEAgPPudkZOvzDJ1yX5i6WUN81/tpK8IMmXllJ2k3zp/O/UWt+W5GVJ3p7k55N8U631I6fa/SkYDAbp9XoL1/V6vWxsbKy4o9UzegwAALTtdmbr3cni40iT5EtucJvnJXneCfpameFwmKtXry5ct7a2luFwuOKOVs/oMQAA0LbbOub0POr3+9ne3k6/3z8cQe31eofLL8JkSEaPAQCAtpVa7zOBbmsuX75cr1271sp97+/vZzQaZTweZ2NjI8Ph8EIE0+TgmNNLly59zIzF1/X7/Qs1YzEAANCcUsrdtdbLi9ad6Dyn59H6+nquXLnSdhutuD5KfPxcr2traxdm9BgAAGiXcEqSZHNzM3t7exd29BgAAGiXcMqhizx6DAAAtOvCT4gEAABA+4yc0lmTySSj0Si7u7sZDAYZDofp9/tttwUAAJyAcEon7ezs3GcCp6tXr2Z7ezubm5tttwcAANwmu/XSOZPJJFtbW5lMJplOp0mS6XR6uHx/f7/lDgEAgNslnNI5o9Eos9ls4brZbJbRaLTijgAAgGXZrZfO2d3dPRwxPW46nWY8Hq+4I+B2OF4cAFhEOKVzBoNBer3ewoDa6/WysbHRQlfArXC8OABwI6XW2nYPhy5fvlyvXbvWdhuccZPJJJcuXcpkMrnPun6/n729vayvr7fQGXAztl0AoJRyd6318qJ1jjmlc/r9fra3t9Pv99Pr9ZIcjJheX+7DLZxNjhcHAG7Gbr100ubmZvb29jIajTIej7OxsZHhcCiYwhnmeHEA4GaEUzprfX09V65cabsN4BY5XhwAuBm79QKwEsPhMGtri9921tbWMhwOV9wRAHCWCKcArITjxQGAm7FbLwAr43hxAOBGhFMAVsrx4gDAInbrBQAAoHXCKQAAAK0TTgEAAGidcAoAAEDrhFMAAABaJ5wCAADQOuEUAACA1gmnAAAAtE44BQAAoHXCKQAAAK0TTgEAAGidcAoAAEDrhFMAAABaJ5wCAADQOuEUAACA1gmnAAAAtE44BQAAoHXCKQAAAK0TTgEAAGidcAoAAEDrhFMAAABaJ5wCAADQOuEUAACA1t3RdgNwEUwmk4xGo+zu7mYwGGQ4HKbf77fdFgAAnBnCKTRsZ2cnW1tbmc1mmU6n6fV6uXr1ara3t7O5udl2ewAAcCbYrRcaNJlMsrW1lclkkul0miSZTqeHy/f391vuEAAAzgbhFBo0Go0ym80WrpvNZhmNRivuCAAAzibhFBq0u7t7OGJ63HQ6zXg8XnFHAABwNgmn0KDBYJBer7dwXa/Xy8bGxoo7AgCAs0k4hQYNh8OsrS3ezNbW1jIcDlfcEQAAnE3CKTSo3+9ne3s7/X7/cAS11+sdLl9fX2+5QwAAOBucSgYatrm5mb29vYxGo4zH42xsbGQ4HAqmAABwhHAKK7C+vp4rV6603QYAAJxZdusFAACgdcIpAAAArRNOAQAAaJ1wCgAAQOtuOZyWUn6klPLBUspbjyz7zlLK75RS3jT/2Tqy7rmllHEp5V2llCefduMAAACcH7czcvovkzxlwfLvq7U+fv6znSSllMcmeXqSx81v86JSygOWbRYAAIDz6ZbDaa31V5L83i1e/alJXlpr/VCt9d1JxkmeeIL+AAAAuABO45jTby6l/OZ8t98HzZddSvK+I9e5Z74MAAAA7mPZcPpDST49yeOTvD/JP50vLwuuWxcVKKU8s5RyrZRy7d57712yHQAAALrojmVuXGv9wPXLpZR/keRV8z/vSfKoI1d9ZJK9G9R4cZIXJ8nly5cXBlhgsclkktFolN3d3QwGgwyHw/T7/bbbAgCA27ZUOC2lPKLW+v75n1+V5PpMvq9M8hOllO9N8qlJBklev8x9AR9rZ2cnW1tbmc1mmU6n6fV6uXr1ara3t7O5udl2ewAAcFtuOZyWUv5NkicleUgp5Z4k35HkSaWUx+dgl93fTvK3kqTW+rZSysuSvD3Jh5N8U631I6faOVxgk8kkW1tbmUwmh8um02mSZGtrK3t7e1lfX2+rPQAAuG23HE5rrV+zYPEP3+T6z0vyvJM0BdzcaDTKbDZbuG42m2U0GuXKlSsr7goAAE7uNGbrBVZsd3f3cKT0uOl0mvF4vOKOAABgOcIpdNBgMEiv11u4rtfrZWNjY8UdAQDAcoRT6KDhcJi1tcWb79raWobD4Yo7AgCA5Qin0EH9fj/b29vp9/uHI6i9Xu9wucmQAADomqVOJQO0Z3NzM3t7exmNRhmPx9nY2MhwOBRMAQDoJOEUOmx9fd2svAAAnAt26wUAAKB1wikAAACtE04BAABonXAKAABA64RTAAAAWiecAgAA0DrhFAAAgNY5zykcM5lMMhqNsru7m8FgkOFwmH6/33ZbAABwrgmncMTOzk62trYym80ynU7T6/Vy9erVbG9vZ3Nzs+32AADg3BJOYW4ymWRrayuTyeRw2XQ6TZJsbW1lb28v6+vrbbW3UkaPAQBYNeEU5kajUWaz2cJ1s9kso9EoV65cWXFXq2f0GACANpgQCeZ2d3cPR0qPm06nGY/HK+5o9Y6OHl9/LKbT6eHy/f39ljsEAOC8Ek5hbjAYpNfrLVzX6/WysbGx4o5W71ZGjwEAoAnCKcwNh8OsrS3eJNbW1jIcDlfc0eoZPYbum0wmueuuu/LsZz87d91118ccRw8AZ5ljTmGu3+9ne3v7Psdbrq2tZXt7+0JMhnR99HhRQL0oo8fQZY4ZB6DLSq217R4OXb58uV67dq3tNrjg9vf3MxqNMh6Ps7GxkeFweCGCaXIw4nLp0qWFIy39fv9CzVgMXWP7BaALSil311ovL1pn5BSOWV9fvxCz8i5i9Bi6y4zjAHSdcAp8jM3Nzezt7V3Y0WPoKseMA9B1wilwHxd59Bi6yjHjAHSd2XoB4Bww4zgAXSecAsA5cP2Y8X6/f3jO5l6vd7jcrvkAnHV26wU6bzKZZDQaZXd3N4PBIMPhMP1+v+22YOUcMw5AlzmVDNBpi87reH12Yed1BAA4W252Khm79QKdNZlMsrW1lclkcjgJzHQ6PVy+v7/fcocAANwq4RTorFs5ryMAAN0gnAKd5byOAADnh3AKdNb18zou4ryOAADdIpwCneW8jgAA54dTyQCddf38jTeardfpMwAAunPaPaeSATpvf3/feR0BABY4a6fdu9mpZIRTAACAc2gymeTSpUuZTCb3Wdfv97O3t7fyL/Sd5xQAAOCC6dpp94RTAACAc6hrp90TTgEAAM6hrp12zzGnAADnTFdm5gSa1bVjTp1KBgDgHFk0M+fVq1dbm5kTaE/XTrtn5BQA4Jw4i6MkQPvO0mn3jJwCAFwAtzIz55UrV1bcFdC29fX1Tmz7JkQCADgnujYzJ8BRwikAwDnRtZk5AY5yzCkAwDnhmFOOM3MzZ41jTgEALoCuzcxJs8zcTNcYOQUAOGfO0syctMMoOmeVkVMAgAukKzNz0hwzN9NFJkQCAIBzxszNdJFwCgAA54yZm+ki4RQAAM6Z4XCYtbXFH/XX1tYyHA5X3BHcP+EUAADOmeszN/f7/cMR1F6vd7jcZEicRSZEAgCAc2hzczN7e3tmbqYzhFMAADinzNxMl9itFwAAgNbdcjgtpfxIKeWDpZS3Hln24FLKa0opu/PfDzqy7rmllHEp5V2llCefduMAAACcH7czcvovkzzl2LLnJHltrXWQ5LXzv1NKeWySpyd53Pw2LyqlPGDpbgEAADiXbjmc1lp/JcnvHVv81CQvmV9+SZKvPLL8pbXWD9Va351knOSJy7UKAADAebXsMacPr7W+P0nmvx82X34pyfuOXO+e+TIAAAC4j6YmRCoLltWFVyzlmaWUa6WUa/fee29D7QAAAHCWLRtOP1BKeUSSzH9/cL78niSPOnK9RybZW1Sg1vriWuvlWuvlhz70oUu2AwAAQBctG05fmeQZ88vPSPKKI8ufXkp5YCnlMUkGSV6/5H0BAABwTt1xq1cspfybJE9K8pBSyj1JviPJC5K8rJRyJcl7k3x1ktRa31ZKeVmStyf5cJJvqrV+5JR7B+icyWSS0WiU3d3dDAaDDIfD9Pv9ttsCAGhdqXXhoaCtuHz5cr127VrbbQA0YmdnJ1tbW5nNZplOp+n1ellbW8v29nY2Nzfbbg8AoHGllLtrrZcXrWtqQiQAjphMJtna2spkMsl0Ok2STKfTw+X7+/stdwgA0C7hFGAFRqNRZrPZwnWz2Syj0WjFHQEAnC3CKcAK7O7uHo6YHjedTjMej1fcEQDA2SKcAqzAYDBIr9dbuK7X62VjY2PFHQEAnC3CKcAKDIfDrK0tfsldW1vLcDhccUcAAGeLcAqwAv1+P9vb2+n3+4cjqL1e73D5+vp6yx0CALTrls9zCsByNjc3s7e3l9FolPF4nI2NjQyHQ8EUACDCKcBKra+v58qVK223AQBw5titFwAAgNYJpwAAALROOAUAAKB1wikAAACtE04BAABonXAKAABA64RTAAAAWiecAgAA0DrhFAAAgNbd0XYDACxvMplkNBpld3c3g8Egw+Ew/X6/7bYAAG6ZcArQcTs7O9na2spsNst0Ok2v18vVq1ezvb2dzc3NttsDALgldusF6LDJZJKtra1MJpNMp9MkyXQ6PVy+v7/fcocAALdGOAXosNFolNlstnDdbDbLaDRacUcAACcjnAJ02O7u7uGI6XHT6TTj8XjFHQEAnIxwCtBhg8EgvV5v4bper5eNjY0VdwQAcDLCKUCHDYfDrK0tfilfW1vLcDhccUcAACcjnAJ0WL/fz/b2dvr9/uEIaq/XO1y+vr7ecocAALfGqWQAOm5zczN7e3sZjUYZj8fZ2NjIcDgUTAGAThFOAc6B9fX1XLlype02AABOzG69AAAAtE44BQAAoHXCKQAAAK0TTgEAAGidcAoAAEDrhFMAAABaJ5wCAADQOuc5BVhgMplkNBpld3c3g8Egw+Ew/X6/7bYAAM4t4RTgmJ2dnWxtbWU2m2U6nabX6+Xq1avZ3t7O5uZm2+2tlJAOAKxKqbW23cOhy5cv12vXrrXdBnCBTSaTXLp0KZPJ5D7r+v1+9vb2sr6+3kJnq7copK+trV3IkA4AnI5Syt211suL1jnmFOCI0WiU2Wy2cN1sNstoNFpxR+2YTCbZ2trKZDLJdDpNkkyn08Pl+/v7LXcIAJw3winAEbu7u4dh7LjpdJrxeLzijtohpAMAqyacAhwxGAzS6/UWruv1etnY2FhxR+0Q0gGAVRNOAY4YDodZW1v80ri2tpbhcLjijtohpAMAqyacAhzR7/ezvb2dfr9/GM56vd7h8osyGZKQDgCsmlPJAByzubmZvb29jEajjMfjbGxsZDgcXphgmnw0pN9ott6L9FgAAKvhVDIA3ND+/v6FDukAwOm62alkjJwCcEPr6+u5cuVK220AABeAY04BAABonXAKAABA64RTAAAAWiecAgAA0DrhFAAAgNaZrReAc2EymWQ0GmV3dzeDwSDD4TD9fr/ttgCAWyScAtB5Ozs72draymw2y3Q6Ta/Xy9WrV7O9vZ3Nzc222wMAboHdegHotMlkkq2trUwmk0yn0yTJdDo9XL6/v99yhwDArRBOAei00WiU2Wy2cN1sNstoNFpxRwDASQinAHTa7u7u4YjpcdPpNOPxeMUdAQAnIZwC0GmDwSC9Xm/hul6vl42NjRV3BACchHAKQKcNh8OsrS1+O1tbW8twOFxxRwDASZxKOC2l/HYp5S2llDeVUq7Nlz24lPKaUsru/PeDTuO+AOCofr+f7e3t9Pv9wxHUXq93uHx9fb3lDgGAW3Gap5L54lrr7x75+zlJXltrfUEp5Tnzv599ivcHAEmSzc3N7O3tZTQaZTweZ2NjI8PhUDAFgA5p8jynT03ypPnllyR5XYRTABqyvr6eK1eutN0GAHBCp3XMaU3y6lLK3aWUZ86XPbzW+v4kmf9+2KIbllKeWUq5Vkq5du+9955SOwAAAHTJaY2cfmGtda+U8rAkrymlvPNWb1hrfXGSFyfJ5cuX6yn1AwAAQIecyshprXVv/vuDSX4myROTfKCU8ogkmf/+4GncFwAAAOfP0uG0lNIrpfSvX07yZUnemuSVSZ4xv9ozkrxi2fsCAADgfDqN3XofnuRnSinX6/1ErfXnSylvSPKyUsqVJO9N8tWncF8AAACcQ0uH01rrbyX5nAXL/2uSL1m2PgAAAOffac3WCwAAACcmnAIAANA64RQAAIDWCacAAAC0TjgFAACgdcIpAAAArRNOAQAAaJ1wCgAAQOuEUwAAAFonnAIAANA64RQAAIDWCacAAAC0TjgFAACgdcIpAAAArRNOAQAAaJ1wCgAAQOuEUwAAAFonnAIAANA64RQAAIDWCacAAAC07o62GwCAs2wymWQ0GmV3dzeDwSDD4TD9fr/ttgDg3BFOAeAGdnZ2srW1ldlslul0ml6vl6tXr2Z7ezubm5tttwcA54rdegFggclkkq2trUwmk0yn0yTJdDo9XL6/v99yhwBwvginALDAaDTKbDZbuG42m2U0Gq24IwA434RTAFhgd3f3cMT0uOl0mvF4vOKOAOB8E04BYIHBYJBer7dwXa/Xy8bGxoo7AoDzTTgFgAWGw2HW1ha/Ta6trWU4HK64IwA434RTAFig3+9ne3s7/X7/cAS11+sdLl9fX2+5QwA4X5xKBgBuYHNzM3t7exmNRhmPx9nY2MhwOBRMAaABwikA3MT6+nquXLnSdhsAcO7ZrRcAAIDWCacAAAC0TjgFAACgdcIpAAAArRNOAQAAaJ1wCgAAQOuEUwAAAFonnAIAANA64RQAAIDW3dF2AwBwEU0mk4xGo+zu7mYwGGQ4HKbf77fdFgC0RjgFgBXb2dnJ1tZWZrNZptNper1erl69mu3t7WxubrbdHgC0wm69ALBCk8kkW1tbmUwmmU6nSZLpdHq4fH9/v+UOAaAdwikArNBoNMpsNlu4bjabZTQarbgjADgbhFMAWKHd3d3DEdPjptNpxuPxijsCgLNBOAWAFRoMBun1egvX9Xq9bGxsrLgjADgbhFMAWKHhcJi1tcVvv2traxkOhyvuCADOBuEUAFao3+9ne3s7/X7/cAS11+sdLl9fX2+5QwBoh1PJAMCKbW5uZm9vL6PRKOPxOBsbGxkOh4IpABeacAoALVhfX8+VK1fabuOWTCaTjEaj7O7uZjAYZDgcpt/vt90WAOeMcAoA3NDOzk62trYym80ynU7T6/Vy9erVbG9vZ3Nzs+32ADhHHHMKACw0mUyytbWVyWRyePqb6XR6uHx/f7/lDgE4T4RTAGCh0WiU2Wy2cN1sNstoNFpxRwCcZ8IpALDQ7u7u4YjpcdPpNOPxeMUdAXCeCacAwEKDweDwdDfH9Xq9bGxsrLgjAM4z4RQAWGg4HGZtbfFHhbW1tQyHwxV3BMB5JpwCAAv1+/1sb2+n3+8fjqD2er3D5c7LCsBpavxUMqWUpyT5/iQPSHJXrfUFTd8nAHA6Njc3s7e3l9FolPF4nI2NjQyHQ8EUgFPXaDgtpTwgyQ8m+dIk9yR5QynllbXWtzd5vwDA6VlfX8+VK1fabgOAc67p3XqfmGRca/2tWusfJ3lpkqc2fJ8AAAB0TNO79V5K8r4jf9+T5PNveO13vSt50pMabgkAAICzpumR07JgWf2YK5TyzFLKtVLKtT/5kz9puB0AAADOoqZHTu9J8qgjfz8yyd7RK9RaX5zkxUly+fLlmte9ruGWAAAAaEVZNH55oOmR0zckGZRSHlNK+fgkT0/yyobvEwAAgI5pdOS01vrhUso3J/mFHJxK5kdqrW9r8j4BgIttMplkNBpld3c3g8Egw+Ew/X6/7bYAuB+l1nr/11qRy5cv12vXrrXdBgDQUTs7O9na2spsNst0Ok2v18va2lq2t7ezubnZdnsAF14p5e5a6+VF65rerRcAYCUmk0m2trYymUwynU6TJNPp9HD5/v5+yx0CcDPCKQBwLoxGo8xms4XrZrNZRqPRijsC4HYIpwDAubC7u3s4YnrcdDrNeDxecUcA3I6mTyUDALASg8EgvV5vYUDt9XrZ2NhooavzxWRTQJNMiAQAnAuTySSXLl3KZDK5z7p+v5+9vb2sr6+30Nn5YLIp4DSYEAkAOPf6/X62t7fT7/fT6/WSHIyYXl8umJ6cyaaAVbBbLwBwbmxubmZvby+j0Sjj8TgbGxsZDoenEkwv8i6ttzLZ1JUrV1bcFW26yNsDzRFOAYBzZX19/dSD0qJdWq9evXphdmk12RRHXfTtgebYrRcA4Cbs0vrRyaYWMdnUxWJ7oEnCKQDATTh/ajIcDrO2tvhj49raWobD4Yo7oi22B5oknAIA3IRdWk02xUfZHmiSY04BAG7C+VMPNDnZFN3Rxe3B5E3d4TynAAA34fyp8FFd2x6cn/fsudl5To2cAgDcxPVdV2/0AXeZD+JGdJrjsW1Gk9vDaTs6edN110d8t7a2zlyQxsgpAMAt2d/fP9VdWo3oNMdj27zT3h6acNddd+VZz3rWDXdB/v7v//6lTjvlC5CTudnIqXAKALBiXds1sks8tlz37Gc/Oy984QtvuP45z3lOnv/855+oti9ATu5m4dRsvQAAK+Z0HM3x2HJdU+fnda7X5ginAAAr5nQczfHYcl1T5+dt+guQyWSSu+66K89+9rNz1113LdwL4LwyIRIAwIp18XQcXeGx5bqmJm9q8guQRbsLX7169cLsLuyYUwCAFXNcZHOafmxNgtM9pz15U1MTLV2U1wUTIgEAnDEmVGlOU4+t/zOS5kJk07MLnxXOcwoAcMZsbm5mb2/vzJ+Oo4uaeGydM5Pruri7cFcIpwAALVlfXz8XIyFn0Wk/trcyCY7/y4ujiS9AHC8tnAIAwP0yqsVxp/0FyHA4zNWrVxeuW2Z24S4RTgEA4H40PaploiWa2l24S0yIBAAA96PJmVRNtMRRpz278Fljtl4AAFhSEyHyopw+BK4zWy8AACypiUlwTLQEHyWcAgDALTrtSXBMtAQftdZ2AwAAcFFdn2hpkYty+hC4TjgFAICWDIfDrK0t/kh+UU4fAtcJpwAA0JLrpw/p9/uHI6i9Xu9wucmQuEgccwoAAC1qYqIl6CLhFAAAWnbaEy1BF9mtFwAAgNYJpwAAALROOAUAAKB1wikAAACtE04BAABonXAKAABA64RTAAAAWiecAgAA0DrhFAAAgNYJpwAAALROOAUAAKB1d7TdAABwMU0mk4xGo+zu7mYwGGQ4HKbf77fdFgC36LRfx0ut9RTbW87ly5frtWvX2m4DAGjYzs5Otra2MpvNMp1O0+v1sra2lu3t7WxubrbdHgD346Sv46WUu2utlxeuE04BgFWaTCa5dOlSJpPJfdb1+/3s7e1lfX29hc4AuBXLvI7fLJw65hQAWKnRaJTZbLZw3Ww2y2g0WnFHANyOpl7HhVMAYKV2d3cznU4XrptOpxmPxyvuCIDb0dTruHAKAKzUYDBIr9dbuK7X62VjY2PFHQFwO5p6HXfMKQCwUo45Beg2x5wCAOdCv9/P9vZ2+v3+4TfvvV7vcLlgCnC2NfU6buQUAGjF/v5+RqNRxuNxNjY2MhwOBVOADjnJ67hTyQAAANA6u/UCAABwpgmnAAAAtG6pcFpK+c5Syu+UUt40/9k6su65pZRxKeVdpZQnL98qAAAA59Udp1Dj+2qt/+ToglLKY5M8Pcnjknxqkl8spXxGrfUjp3B/AAAAnDNN7db71CQvrbV+qNb67iTjJE9s6L4AAADouNMIp99cSvnNUsqPlFIeNF92Kcn7jlznnvkyAAAAuI/7DaellF8spbx1wc9Tk/xQkk9P8vgk70/yT6/fbEGpheesKaU8s5RyrZRy7d577z3ZvwIAAIBOu99jTmutf+lWCpVS/kWSV83/vCfJo46sfmSSvRvUf3GSFycH5zm9lfsCAOD8mEwmGY1G2d3dzWAwyHA4TL/fb7stVszzgFLryfNgKeURtdb3zy9/a5LPr7U+vZTyuCQ/kYPjTD81yWuTDO5vQqTLly/Xa9eunbgfAAC6ZWdnJ1tbW5nNZplOp+n1ellbW8v29nY2Nzfbbo8V8Ty4OEopd9daLy9ct2Q4/fEc7NJbk/x2kr91JKx+e5JvSPLhJM+qtf7c/dUTTgEALo7JZJJLly5lMpncZ12/38/e3l7W19db6IxV8jy4WG4WTpeaEKnW+nW11v+h1vrZtdavuB5M5+ueV2v99Frrn7uVYAoAwMUyGo0ym80WrpvNZhmNRivuiDZ4HnBdU6eSAQCAm9rd3c10Ol24bjqdZjwer7gj2uB5wHXCKQAArRgMBun1egvX9Xq9bGxsrLgj2uB5wHVLHXN62hxzCgBwcTjWkMTz4KJp7JhTAAA4qX6/n+3t7fT7/cORs16vd7hcILkYPA+4zsgpAACt2t/fz2g0yng8zsbGRobDoUByAXkeXAyNnUrmtAmnAAAA55fdegEAADjThFMAAABaJ5wCAADQOuEUAACA1gmnAAAAtE44BQAAoHXCKQAAAK0TTgEAAGidcAoAAEDrhFMAAABaJ5wCAADQOuEUAACA1gmnAAAAtE44BQAAoHXCKQAAAK0TTgEAAGidcAoAAEDrhFMAAABaJ5wCAADQOuEUAACA1pVaa9s9HCql3JvkPW33AQAAQCM+rdb60EUrzlQ4BQAA4GKyWy8AAACtE04BAABonXAKAABA64RTAAAAWiecAgAA0DrhFAAAgNYJpwAAALSuM+G0lLJWSlmbX/74UsoTSikPPqXaH7dg2UNOo/aCup+5xG1LKeXzSyl/pZTyVfPL5ZT6ujyv+eXL9His5qNLKZ88v3xnKeVppZTPWrLmZ59GbwvqfvzRx7KU8sWllG8rpfzlU7yP9fnz9pNPq+aR2n/7lOp0fjtb9vnb5HbGxyql/Jm2ezgvVv28PY33iSZeE5p8DTt2P6fymnuk3qm8PzT1XtbUe++89ql/VpjXauJzzZNLKT9USnllKeUV88tPOY3aN7i/v7/k7Z9cSrlSSrnz2PJvWKJmKaX8tVLKV88vf0kp5Z+XUv729W3vNJRSfukUajzk2N//+7zXZy7z+jh/Xj14fvmhpZQfK6W8pZQyKqU8com631tK+cKT3v4GNR9cSvn7pZRvnP9/fXsp5VWllH9cSnnQkrWb2R5qrWf+J8lXJvlAkvcneWqSX0/yS0nuSfLlS9T94nmNe5O8OsmdR9a9saF/y3tPeLsvSzJO8nNJ7pr//Px82Zct0c8XJbmW5BeT/H6SVyX51SSvS/KoJeo+J8m7k7wzyTfOf/9wkrclubpE3Y/M/83fk+Sxp/j/8uYkD5pf/r+T/IckfzfJa5I8/4Q1X3Tk8maS9yb55STvS7K1RK9Xj/18W5Lfvf73EnXPxXZ20m1sftumtrP1JN89f/7/t/lj8R+TfP2S/9bL8+fUv0ryqPnz9b8leUOSz12i7hvnz/9PP8X/lxckeciRvn9r/ri+J8kXLVH3KUcuf9L8deY3k/xEkocvUfdTkvxQkh9M8meSfGeStyR5WZJHLFH3E5M8P8mPJ/naY+tedMKajTxv7+c+l9nOGnlNaPA17NRfc9Pc+8Opv5fNazX13nvqnxXS3Oeaf5ZkO8nT5/9nm/PL20m+/7Qek2P3ucx29g+T/Mq87/+c5FuOrFtmO3tRkp9K8socvPf8ZJK/nuSlJ30ccvCaffTnLUk+dP3vJXp945HLfzfJLyR5xrzn71ui7tuPXB4l+dYkj0zy9Ules0Tde+fP3fckeWGWeB8/UnM7yT/KwfvZ65L8QJK/kIPPJK9Yom5j28NS/+BV/ST5jRx8UHhMkj9M8ufmyz8tybUl6r4hyePml5+WZDfJF1y/zyXq/vMb/PxAkj88Yc135Mgb+JHlj0nyjiUf24ceqfUz88tfmuTVS9R9W5JPyMGHusmR++gleeuS/X5Wkufl4I3yzTl4c7vPY3Obdd965PK1JJ8wv3zHSV8Yj70o/nKSJ8wv/9kln7eT+Yvh30/yHfOf379+ecnHthPbWRPb2LxuU9vZK3LwpvXIHHyg/XtJBklekuQfLlH39Un+cpKvycGH2qfNl39Jkl9bou67k/yTHHxgfn0O3ng/9aT15jXfcuTyLyf5vPnlz1jy+XV0O7sryT+YP2e/NcnPLlH355N8y/z15TeTPDvJo+fLlnlDf3kOgvpX5uDD3cuTPPD4v+WMPG+b2s6aeu9t6jXs1F9z09z7w6m/lx15bJt47z31zwpp7nPNf7rB8pJkd4m6f3iDn0mSDy9R9y1J7phf/uQchIbvu/4YLVN3/vvjkvzXJB9/5Dn2lhPWvB50P3O+vd6Zg/e0T0vyacs8b49cfmOS3pHeT9Tr/PbvOnL57mPr3rRsvzn4fPD35tvHO+evNZ9xwppvqh99nv7OKfbayPZQa+3Obr211v9Sa313Dr5Fetd82Xuy3K7JH19rfdu81k/l4MPCS0opX5WkLlH3byR5a5K7j/1cS/LHJ6x5Rw6+/T3ud3KwkZ3UA2qt984vvzcHLwSptb4myaUl6n6k1vrfk/xBkv+egxew1FqnS9Scl6hvrbV+e611I8nfTPKwJP++lPIflqj7h0d2I/rdJH9qfvmOnM7u759Ya31jktRafyvJA5ao9bj57XtJ/nGt9buS/H6t9bvml0+sQ9tZE9tY0tx2dmet9V/WWu+ptX5vkq+ote7m4N/xV5ao+3G11p+rtf6bHGwbP5WDC6/NR5/DJ/H7tdb/q9b66ByMEg2SvLGU8sullGeetNdSyh3zy59Qa33DvNf/lOSBS/R61OVa69+ttb6n1vp9OfiAc1IPr7X+QK31BUk+udb6j2qt7621/kDmr5Mn9Om11ufUWn+21voVOfjA9EtL7uLc1PO2qe2sqffepl7DGnvNnTvN94em3suaeu9t4rNCU59r/qiU8sQFyz8vyR8tUfcPkgxqrZ947Kefg70ATuqOWuuHk6TW+gdJvjzJJ5ZSfjLJxy9R93rNP0nyhlrrH8///nAORthv2/y18OVJXpzkc2qtv53kT+av5e9ZotdPKKV8binlf8zB82J6pPcT9Tr3ulLKd5dSPmF++SuTg93oc7D30knVeX+7tdbvqbU+Lslfy8F2vH3Cmmvz3XcflWT9+i7e8/ecZZ4HTW0PueP+r3I2lFLWaq2zJN9wZNkDstwD+yellE+ptf6XJKm1vq2U8iU52AXk05eo+4YcfON3nxfsUsp3nrDmjyR5QynlpTn4Nik5eKI9PQe7wJzUtVLKDyd5bQ52g3rdvM8/neXeIN9YSvmJHLyZvzYHHzx+PslfTPL2Jep+zDECtdbXJ3l9KeXbkvwvS9T9P5L861LKm5N8MAePy79L8tk52DXmJD6zlPKb857vLKU8qNb6+/NjMk78gbHW+t4kTyulPDXJa0op33fSWsd1aDtrYhtLmtvOpqWUzVrrTinly5P8XpLUWmfLHPeSgzeHL8vB7qy1lPKVtdafLaV8UZZ74z1Ua/33OfgA+i05GHkY5uADxO36wSTbpZQXJPn5Uso/S/LTORjlfdMSLT6slHI1B9vZJ5ZSSp1/fZvlPowfve2P3WTd7Xrgke0stdbnlVLuycHud+snrNnU87ap7ayp995GXsMaes1t5P0hzbyXJc299zbxWaGpzzVfn+SHSin9fPTLoEflYJTz65eo+2M5CNAfWLDuJ5ao+59LKV9Ua/13SVJr/UiSK6WUf5Dkry5R97+UUtZrrfu11sPjC0spn5IlvrSqtf5MKeXVSb6nlPKNWe5zx3XvT/K988u/V0p5RK31/fNg9uEl6n5zkm9P8q75399aSpkm+bdJvm6Juvf5PFBrvb6r83NPWPP5ORh9TQ5eF+8qpdQkj02yzJdrX59mtoeUj76Hn12llM/LwfD7Hx1bfmeSzVrrvzph3b+U5N5a65uPLf+kJN9ca33eCes+OMkf1Vr/v5Pc/iZ1H5vkK3LwzV/JwZPhlbXWE4e9cjAhxd/MwZP0zUl+pNb6kfm3QQ876TdW8xGSr87Bt0A/leSJSb42B99i/uBJvxUtpXxtrXWZF+ub1X5ADo7d+ox8dCTiF+bfOJ6k3vHRlb1a65+UgwP0/5da608v0+/8Pv50Dl5cPr/WuswHhE5tZ01tY/Pafz4HH2hOczv77BzscvoZORiJ+oZa638qpTw0ydfUWv/5Cet+Tg6OS5nlYDfW/zMHx9P8TpK/uShU3GLdl9Zan36S295P3SfloMej29jPJPnR+TfZJ6n5HccWvajWeu/8g9ILa61//YR1v3t++/1jyzeSvKDW+rQT1n1hDnYt/MVjy5+S5AdqrYMT1m3iedvUe1lT772NvIYdq3Uqr7lNvj+c9nvZvGYj771NfFZo6nPNkfqfkiPb2fUvWc6a+b8385Hp4+su1Vp/55Tvr5eD3WY/eAq1PifJ/1Rr/X+X72xh/Qfk4HCKpV/b5q9bd9Ra/+sp1Fo//p5zGub/3lJr/fB8m3t8DnbxXWZk/nrtU98eOhFOAYDllFJeXmtdZsRE3RXX7VKvXat7mjVLKZ9Za33n/V9T3bNQU93Tq1lK+bjjXy6XUh5Sa/3dk9bsxDGnpZRPLKU8v5Ty46WUrz227kUXoe6xml9zWr3ez33+nLrN1O1Sr12r25VeSymvOs16XazbpV67WHeBP6tu5+p2qdeu1T3Nmq8+xVrqNl9T3SVrloPTVN2TZK+U8urysacrWqrXrhxz+qM5mM3v5Um+oZTyV3MwBf+HknzBBal7vObTTqPXUsoTbrQqB8P+6p6wbpd67VrdLvV6E8tMzHFe6nap1y7WPa6pXaXUba5ul3rtWt3bqllKudHhFyUHs+GeiLrd6rVrdZvqNQeHFD15PmfA03JwPP7X1Vr/47z2iXUlnH76kV0vfraU8u05mN3wKy5Q3aZ6fUOSf5fFT6RPVnepul3qtWt1u9TrjfzGKdfrYt0u9drFusDp+Rs5mL38QwvWfc2CZeq2W1Pd5momx2ZdL6W8I8lPl1Kek2W/TKpLnIdmVT85OIfb2rFlz8jB+X/ecxHqNtjrW3Mwhfmide9T9+R1u9Rr1+p2qdfbvP+XX/S6Xeq1a3WzxLkN1W2nbpd67Vrd262Z5JeS/M83WPfuJfq48HW71GvX6jbY67Ukn3Js2SNzMPv+5KR1a62dGTn9tzmYVvxwdsNa60tKKR/IwcnAL0Ldpnr9ztz42ONvUXepuk3UVLe5mk3WvVVdOl6rqbpd6rVrdZ/dQE11m63bpV67Vvd2az4tNzh/Y631MUv0oW63eu1a3aZ6fU6Shyc5nJ231npPOTiV3TcvUddsvQBnRSnljbXWGx33eiHqdqnXs1K3lPKW3GQ3qlrrZ5+wB3UbqtulXrtWt6leb+P+OzNjcdfqdqnXrtU9S712ZeT0Pkopr6q1/q8XuW6XelW3uZrqNlezybpwiq4/P79p/vvH57//tyTLnMdP3ebqdqnXrtVtqtdb1aW9KbpWt0u9dq3umem1s+E03Zs10UyX6jZVU93majZZd5GlZrg7J3W71OuZqFtrfU+SlFK+sNb6hUdWPaeU8qtJvvskDajbXN0u9dq1uk31ejstqNtY3S712rW6Z6bXTpzn9Aa6NmuimS7Vbaqmus3VbLLuIl06Xqupul3q9azV7ZVSNq//UUr5wiS9U+hF3ebqdqnXrtVtqlegQefqmNOztL90W3W71Ku6zdVUt7maJ6nbpeO1mqrbpV67WHde+wk5OCf2J83v478l+Ru11qW+YFG3ubpd6rVrdZvq9Rbu9zdqrZ+r7unX7VKvXat7lnrt8m69i5yZ/aVbrNulXtVtrqa6zdU8Sd0uHa/VVN0u9dqpuqWUq0f+/LEcjA5N539/cU44+q9uc3W71GvX6jbV6204S3tTnLe6Xeq1a3XPTK/nLZyemf2lW6zbpV7Vba6mus3VvO26XTpeq6m6Xeq1g3X7899/LsnnJXlFDkaLvjzJr5ykT3Ubr9ulXrtWt5Feb3Wvh1rrq9W9vbpd6rVrdbvU63XnLZwCnGW9UspmrXUnOf3jtTpSt0u9dqJurfW75jVeneQJtdbJ/O/vTPKTJ21Q3ebqdqnXrtVtqtd0aG+KDtbtUq9dq9ulXg/UWs/NT5LfuOh1u9Srut3rtWt1z1qvSZ6Q5M1JfjvJu5O8KcnnnkI/nanbpV67VjfJO5M88MjfD0zyzlPoVd2G6nap167VbbDXX72VZeqejZrqdq/X8zZyemb2l26xbpd6Vbe5muo2V/O263bpeK2m6nap1y7WnfvxJK8vpfxMDna3+qokL1minrrN1+1Sr12r21SvZ35vig7X7VKvXavbmV47MVtv12ZNbKJul3pVt7ma6jZXs+G63zG/ePQYqJL5MVC11m8873W71GsX6x6p/4Qkf2H+56/UU5qZVN3m6nap167VbbDmj6YDMxZ3rW6Xeu1a3U712pFw+mnziwv3a661nmhiii7V7VKv6nav167V7VKvx+q/OslfrR89Bqqf5CdrrU+5KHW71GsX6wLNOLbXQ8nH7vVQa63fq+7J6nap167V7VKvh8W6EE6vK6X8av3Y2Q0XLjvPdbvUq7rN1VS3uZoN131nks+ptX5o/vcDk7y51vqZF6Vul3rtYl2gGV3bm6JLdbvUa9fqdqnX67p2zGln9pdusG6XelW3uZrqNlezybpdOl6rqbpd6rWLdYEG1A7NWNy1ul3qtWt1u9TrdV0Lp9+Q5EdLKR+zX/MFq9ulXtVtrqa6zdVsrG6t9XmllJ/LR4+BWvp4j67V7VKvXawLNO7RSf74yN9/nOROdU+lbpd67VrdzvTaid16S8f2l26ibpd6Vbd7vXatbpd6BeD8KKV8e5K/luToXg+jWuvz1V2ubpd67VrdLvXalZHT/vz30f2aPynz/ZovSN0u9apuczXVba5mk3UBOAe6tjdFl+p2qdeu1e1Sr50YOb2udGzWxCbqdqlXdbvXa9fqdqlXAABubq3tBm5TZ/aXbrBul3pVt7ma6jZXs8m6AADcQFd2672ua7MmmulS3aZqqttczSbrAgBwA53arTdJSilPyEf3a/6V05rdsEt1u9Srus3VVLe5mk3WBQBgsc6FUwAAAM6frh1zCgAAwDkknAIAANA64RQAAIDWCacAAAC0TjgFAACgdf8/C8iypq1pAjwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1152x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "het_eventstudy = PanelOLS(df.y, df[ind_cols],  \n",
    "                         entity_effects = True, time_effects = True,\n",
    "                         check_rank = False, drop_absorbed = True).fit(cov_type = 'robust')\n",
    "\n",
    "het_err_series = het_eventstudy.params - het_eventstudy.conf_int()['lower']\n",
    "het_coef_df = pd.DataFrame({'coef': het_eventstudy.params.values[1:48],\n",
    "                        'err': het_err_series.values[1:48],\n",
    "                        'varname': het_err_series.index.values[1:48]\n",
    "                       })\n",
    "\n",
    "fig, ax = plt.subplots(figsize = (16, 10))\n",
    "het_coef_df.plot(x = 'varname', y = 'coef', kind = 'bar', \n",
    "             ax = ax, color = 'none', \n",
    "             yerr = 'err', legend = False)\n",
    "ax.set_ylabel('')\n",
    "ax.set_xlabel('')\n",
    "ax.axhline(y = 0, color= 'r', linestyle='-')\n",
    "ax.scatter(x = np.arange(het_coef_df.shape[0]), \n",
    "           marker = 'o', s = 50, \n",
    "           y = het_coef_df['coef'], color = 'black')\n",
    "ax.xaxis.set_ticks_position('none')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "932ad457",
   "metadata": {},
   "source": [
    "##### Event study with constant treatment effects, no anticipatory behavior and parallel trends (Sun and Abraham 2020). Droppping two leads."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 164,
   "id": "72d4e3de",
   "metadata": {},
   "outputs": [],
   "source": [
    "leads_2 = {'dd_-24', 'dd_-4', 'dd_-1'}\n",
    "ind_cols_2 = [ind for ind in dd if ind not in leads_2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 170,
   "id": "13b48f50",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA7AAAAJbCAYAAADdSXg3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA5vElEQVR4nO3de5hsaV0f+u/bM+Ax1eUtcnPDgNotBvMIMi3qYScZoij2o6KGpMAcr5OMl+iTccw5YEhmzEUlmKh4JWTUoIlaGkWItggaCW6PHtkzgIBcukWQsSeCRqW63Coz9Z4/qvamp6cbZlf16u63+/N5nn521VrV3/5NT1et+lattarUWgMAAAAn3dJxDwAAAAAPhAILAABAExRYAAAAmqDAAgAA0AQFFgAAgCYosAAAADThUApsKeWHSynvLqW8cdeyjyqlvLKUsjn79yMP+N6nlVLeWkrZKqU85zDmAQAA4PQ5rHdg/3OSp+1Z9pwkv1JrXU3yK7Pr91FKuSbJ9yf53CSPS/KsUsrjDmkmAAAATpFDKbC11lcn+d97Fj89yYtnl1+c5Av3+dYnJdmqtb691vpXSX5y9n0AAABwH10eA/uwWuvdSTL796H73OZcknftun7XbBkAAADcx7XH/PPLPsvqvjcs5aYkNyVJr9e7/hM/8RO7nAsAAIBjcscdd/xRrfUhe5d3WWD/sJTyiFrr3aWURyR59z63uSvJo3Zdf2SS7f3Caq0vSvKiJFlbW6sXL1487HkBAAA4AUop79xveZe7EL8syZfPLn95kpfuc5vXJFktpXxsKeXBSZ45+z4AAAC4j8P6GJ2fSPIbSR5bSrmrlHJjkucleWopZTPJU2fXU0r5mFLKRpLUWu9J8vVJfinJm5P8VK31TYcxEwAAAKfLoexCXGt91gGrPnOf224nWd91fSPJxmHMAQAAwOnV5S7EAAAAcGgUWAAAAJqgwAIAANAEBRYAAIAmKLAAAAA0QYEFAACgCYfyMTon0Wg0ynA4zObmZlZXVzMYDNLv9w8ld21tLZcuXcqtt956aLk33HBDkuRVr3rVwlmXdfU7AAAAOA6l1nrcM1y1tbW1evHixQPXX7hwIevr65lMJhmPx+n1ellaWsrGxkbOnz8/98+9nDsejzOZTA4tt4tS3NWsAAAAXSul3FFrXbvf8tNWYEejUc6dO5fRaHS/df1+P9vb21leXr7qn9lVbhdFs6tZAQAAjsJBBfbUHQM7HA4zmUz2XTeZTDIcDk9M7mg0yvr6ekaj0ZXs8Xh8ZfnOzs6JmRUAAOC4nboCu7m5mfF4vO+68Xicra2tE5PbVdHs6ncAAABwnE5dgV1dXU2v19t3Xa/Xy8rKyonJ7apodvU7AAAAOE6nrsAOBoMsLe3/n7W0tJTBYHBicrsqml39DgAAAI7TqSuw/X4/Gxsb6ff7V8phr9e7snzekxftzr1cDhfN7apodvU7AAAAOE6n7izEl+3s7GQ4HGZraysrKysZDAaHUtx2dnZy/fXX59KlS7ntttsWzu3y4266+h0AAAB06cx8jE6LDrsUAwAAtOygAnvtcQzDfS0vL+etb33rcY8BAABwop26Y2ABAAA4nRRYAAAAmqDAAgAA0AQFFgAAgCYosAAAADRBgQUAAKAJCiwAAABNUGABAABoggILAABAExRYAAAAmqDAAgAA0AQFFgAAgCYosAAAADRBgQUAAKAJCiwAAABNUGABAABoggILAABAExRYAAAAmqDAAgAA0AQFFgAAgCYosAAAADRBgQUAAKAJCiwAAABNUGABAABoggILAABAExRYAAAAmqDAAgAA0AQFFgAAgCYosAAAADRBgQUAAKAJCiwAAABNUGABAABoggILAABAExRYAAAAmqDAAgAA0AQFFgAAgCYosAAAADRBgQUAAKAJCiwAAABNUGABAABoggILAABAExRYAAAAmtBpgS2lPLaU8rpdX+8tpdy85zY3lFL+bNdtbu1yJgAAANp0bZfhtda3JnlCkpRSrknyB0less9Nf63W+nldzgIAAEDbjnIX4s9M8ru11nce4c8EAADglDjKAvvMJD9xwLrPKKW8vpTyi6WUTzrCmQAAAGjEkRTYUsqDk3xBkp/eZ/WdSR5da318ku9N8nMHZNxUSrlYSrn4nve8p7NZAQAAOJmO6h3Yz01yZ631D/euqLW+t9a6M7u8keRBpZSP3ud2L6q1rtVa1x7ykId0PzEAAAAnylEV2GflgN2HSykPL6WU2eUnzWb64yOaCwAAgEZ0ehbiJCml/LUkT03y1buWfU2S1FpfmOQZSb62lHJPkktJnllrrV3PBQAAQFs6L7C11j9P8tf3LHvhrsvfl+T7up4DAACAth3lWYgBAABgbgosAAAATVBgAQAAaIICCwAAQBMUWAAAAJqgwAIAANAEBRYAAIAmKLAAAAA0QYEFAACgCQosAAAATVBgAQAAaIICCwAAQBMUWAAAAJqgwAIAANAEBRYAAIAmKLAAAAA0QYEFAACgCQosAAAATVBgAQAAaIICCwAAQBMUWAAAAJqgwAIAANAEBRYAAIAmKLAAAAA0QYEFAACgCQosAAAATVBgAQAAaIICCwAAQBMUWE610WiU22+/Pc9+9rNz++23ZzQaHVruYx/72Fx33XXN5HbxewAAgKNUaq3HPcNVW1tbqxcvXjzuMTjhLly4kPX19Uwmk4zH4/R6vSwtLWVjYyPnz59fOHc8HmcymTSTe9i/BwAA6Eop5Y5a69r9liuwnEaj0Sjnzp3b953Gfr+f7e3tLC8vy50zFwAAunRQgbULMafScDjMZDLZd91kMslwOJS7QC4AABwHBZZTaXNzM+PxeN914/E4W1tbchfIBQCA46DAciqtrq6m1+vtu67X62VlZUXuArkAAHAcHAPLqdTaMaWt5QIAQJccA8uZ0u/3s7GxkX6/f+UdyF6vd2X5vKVtd+7S0lJTuYf5ewAAgOPgHVhOtZ2dnQyHw2xtbWVlZSWDweBQStvOzk6uv/76XLp0KbfddlsTuV38HgAAoAs+RgcAAIAm2IUYAACApimwAAAANEGBBQAAoAkKLAAAAE1QYAEAAGiCAgsAAEATFFgAAACaoMACAADQBAUWAACAJiiwAAAANEGBBQAAoAkKLAAAAE1QYAEAAGiCAgsAAEATFFgAAACaoMACAADQBAUWAACAJiiwAAAANEGBBQAAoAkKLAAAAE1QYAEAAGhC5wW2lPKOUsobSimvK6Vc3Gd9KaV8Tyllq5Ty26WUJ3Y9EwAAAO259oh+zlNqrX90wLrPTbI6+/q0JD84+xcAAACuOAm7ED89yY/Wqd9M8hGllEcc91AAAACcLEdRYGuSV5RS7iil3LTP+nNJ3rXr+l2zZQAAAHDFUexC/ORa63Yp5aFJXllKeUut9dW71pd9vqfuXTArvzclyXXXXdfNpAAAAJxYnb8DW2vdnv377iQvSfKkPTe5K8mjdl1/ZJLtfXJeVGtdq7WuPeQhD+lqXAAAAE6oTgtsKaVXSulfvpzks5O8cc/NXpbky2ZnI/70JH9Wa727y7kAAABoT9e7ED8syUtKKZd/1o/XWl9eSvmaJKm1vjDJRpL1JFtJ/jzJV3Y8EwAAAA3qtMDWWt+e5PH7LH/hrss1yT/pcg4AAADad1SfAwsAnHKj0SjD4TCbm5tZXV3NYDBIv98/7rEAOEUUWABgYRcuXMj6+nomk0nG43F6vV5uueWWbGxs5Pz58wtld1GMuyrbSjxAt8p0D962rK2t1YsXLx73GABApqXt3LlzGY1G91vX7/ezvb2d5eXlubL3K8ZLS0sLFePLmePxOJPJ5FAyu5oV4KwqpdxRa13bu7zzj9EBAE634XCYyWSy77rJZJLhcDhX7mg0yvr6ekajUcbjcZJkPB5fWb6zs7NQ5uWZF83salYA7k+BBQAWsrm5eaW07TUej7O1tTVXbhfFuKuy3VUuAPelwAIAC1ldXU2v19t3Xa/Xy8rKyly5XRTjrsp2V7kA3JcCCwAsZDAYZGlp/6cUS0tLGQwGc+V2UYy7Kttd5QJwX07iBAAsrIsTGHVxcqiuTjjV5YmsAM6ig07i5GN0AICFnT9/Ptvb2xkOh9na2srKykoGg8FCpa3f72djY+PAYjxPdheZXeYCcF/egQUATrSdnZ1DLcZdZXaZC3DWHPQOrAILAADAieJzYAEAAGiaAgsAAEATFFgAAACaoMACAADQBAUWAACAJiiwAAAANEGBBQAAoAkKLAAAAE1QYAEAAGiCAgsAAEATFFgAAACaoMACAADQBAUWAACAJiiwAAAANEGBBQAAoAkKLAAAAE1QYAEAAGiCAgsAAEATFFgAAACaoMACAADQBAUWAACAJiiwAAAANEGBBQAAoAkKLAAAAE1QYAEAAGiCAgsAAEATFFgAAACaoMACAADQBAUWAACAJiiwAAAANEGBBQAAoAkKLAAAAE1QYAEAAGiCAgsAAEATrj3uAQCOwmg0ynA4zObmZlZXVzMYDNLv9497LAAAroICC5x6Fy5cyPr6eiaTScbjcXq9Xm655ZZsbGzk/PnzC2V3VYy7yFXigaPi8QboSqm1HvcMV21tba1evHjxuMcAGjAajXLu3LmMRqP7rev3+9ne3s7y8vJc2ZeL8Xg8zmQySa/Xy9LS0sLFuIvc/Ur8YcwKsFdrjzctvRAJZ0kp5Y5a69r9liuwwGl2++235+abb854PL7ful6vlxe84AW58cYbrzq3q2LcRW6XJR5gty4fb7oohF2V7dZKPJxEBxVYJ3ECTrXNzc19y2uSjMfjbG1tzZU7HA4zmUz2XTeZTDIcDk9MblezAuzV1ePNhQsXcu7cudx88815/vOfn5tvvjnnzp3LhQsX5p51NBplfX09o9HoynZiPB5fWb6zs3OicoEpBRY41VZXV9Pr9fZd1+v1srKyMlduV8W4i9yuZgXYq4vHm64KYUsvRALvp8ACp9pgMMjS0v4PdUtLSxkMBnPldlWMu8jtalaAvbp4vOmqELb0QiTwfgoscKr1+/1sbGyk3+9feVLV6/WuLJ/3WKyuinEXuV3NCrBXF483XRXCll6IBN5PgQVOvfPnz2d7ezsveMEL8pznPCcveMELsr29vdCJNHYX48tP1g6jGHeR21WJB9iri8ebrgphSy9EAu/nLMQAC9jZ2clwOMzW1lZWVlYyGAwOpRB2kdvVrAB7HebjzVF8HJqzEMPJ42N0AABoUpeFsKUXIuEsUWABAGiWQghny0EF9trjGAYAAK7G8vJybrzxxuMeAzhmTuIEAABAEzotsKWUR5VSfrWU8uZSyptKKf90n9vcUEr5s1LK62Zft3Y5EwAAAG3qehfie5J8U631zlJKP8kdpZRX1lp/Z8/tfq3W+nkdzwIAAEDDOn0HttZ6d631ztnlUZI3JznX5c8EAADgdDqyY2BLKY9J8ilJ/r99Vn9GKeX1pZRfLKV80lHNBAAAQDuO5CzEpZTlJD+T5OZa63v3rL4zyaNrrTullPUkP5dkdZ+Mm5LclCTXXXddtwMDAABw4nT+Dmwp5UGZltf/Wmv92b3ra63vrbXuzC5vJHlQKeWj97ndi2qta7XWtYc85CFdjw0AAMAJ0/VZiEuSH0ry5lrrdx5wm4fPbpdSypNmM/1xl3MBAADQnq53IX5yki9N8oZSyutmy/55kuuSpNb6wiTPSPK1pZR7klxK8sxaa+14LgAAABrTaYGttV5IUj7Ibb4vyfd1OQcAAADtO7KzEAMAAMAiFFgAAACaoMACAADQBAUWAACAJiiwAAAANEGBBQAAoAkKLAAAAE1QYAEAAGiCAgsAAEATFFgAAACaoMACAADQBAUWAACAJiiwAAAANEGBBQAAoAkKLAAAAE1QYAEAAGiCAgsAAEATFFgAAACaoMACAADQBAUWAACAJiiwAAAANOHa4x4AAAA4fUajUYbDYTY3N7O6uprBYJB+v3/cYx2oq3lb+z2cdKXWetwzXLW1tbV68eLF4x4DAADYx4ULF7K+vp7xeJzJZJJer5elpaVsbGzk/PnzC2V3UQgvzzuZTDIejw9t3q5yz4JSyh211rX7LVdgAQCAwzIajXLu3LmMRqP7rev3+9ne3s7y8vJc2V0Uwq7m7fL3cBYcVGAdAwsAABya4XCYyWSy77rJZJLhcDhX7mg0yvr6ekajUcbjcZJkPB5fWb6zs3Oi5u0q96xTYAEAgEOzubl5pWDuNR6Ps7W1NVduV4Wwq3m7yj3rFFgAAODQrK6uptfr7buu1+tlZWVlrtyuCmFX83aVe9YpsAAAwKEZDAZZWtq/ZiwtLWUwGMyV21Uh7GrernLPOgUWAAA4NP1+PxsbG+n3+1cKZ6/Xu7J83hMXdVUIu5q3q9yzzlmIAQCAQ7ezs5PhcJitra2srKxkMBgsXNq6/FiaLubtMve08zE6AABA8xTCs+GgAnvtcQwDAAAwj+Xl5dx4443HPQbHxDGwAAAANEGBBQAAoAkKLAAAAE1QYAEAAGiCAgsAAEATFFgAAACaoMACAADQBAUWAACAJiiwAAAANEGBBQAAoAkKLAAAAE1QYAEAAGiCAgsAAEATFFgAAACaoMACAADQBAUWAACAJiiwAAAANEGBBQAAoAkKLAAAAE1QYAEAAGiCAgsAAEATFFgAAACaoMACAADQBAUWAACAJiiwAAAANEGBBQAAoAkKLAAAAE24tusfUEp5WpIXJLkmye211uftWV9m69eT/HmSr6i13tn1XACcTaPRKMPhMJubm1ldXc1gMEi/3z/usQDg1DrMbW+ptR7yeLvCS7kmyduSPDXJXUlek+RZtdbf2XWb9STfkGmB/bQkL6i1ftoHyl1bW6sXL17sbG4ATqcLFy5kfX094/E4k8kkvV4vS0tL2djYyPnz5xfKHo1GWVtby6VLl3LrrbcqxgCQ+be9pZQ7aq1r91vecYH9jCTfUmv9nNn1b06SWuu377rNf0zyqlrrT8yuvzXJDbXWuw/KVWABuFqj0Sjnzp3LaDS637p+v5/t7e0sLy/Pld1lMQaAVi2y7T2owHa9C/G5JO/adf2uTN9l/WC3OZfkwAKbt741ueGGw5kQgDNh5+678wvjce7dZ90143FG11+f5Uc84qpz77n33tTf+I287N5dyeNxkqTecEPu/YzPyDXXXDPn1ADQri62vV2fxKnss2zvW74P5DYppdxUSrlYSrn4vve971CGA+Ds+PNLl3LvZLLvunsnk1y6dGmu3Pe8+93JQXsz1Zp3v/vdc+UCQOu62PZ2/Q7sXUketev6I5Nsz3Gb1FpflORFyXQX4rzqVYc6KACn26/efntuvvnmjGfvju7W6/Xygttuy8fdeONV5373s5+d5z//+fuvnEzynC/+4nz7t3/7/usB4BRbaNtb9nufs/t3YF+TZLWU8rGllAcneWaSl+25zcuSfFmZ+vQkf/aBjn8FgHkMBoMsLe2/2VtaWspgMJgrd3V1Nb1eb991vV4vKysrc+UCQOu62PZ2WmBrrfck+fokv5TkzUl+qtb6plLK15RSvmZ2s40kb0+yleQ/Jfm6LmcC4Gzq9/vZ2NhIv9+/sjHt9XpXls97AqeuijEAtG73tvfyi72Lbns7PQtxV5yFGIB57ezsZDgcZmtrKysrKxkMBnOX18sun4V4MplkPB47CzEA7DLPtvdYPkanKwosACdNF8UYAM6q4/oYHQA4E5aXl3PjHCeBAgAeuK5P4gQAAACHQoEFAACgCQosAAAATVBgAQAAaIICCwAAQBMUWAAAAJqgwAIAANAEBRYAAIAmKLAAAAA0QYEFAACgCQosAAAATVBgAQAAaIICCwAAQBMUWAAAAJqgwAIAANAEBRYAAIAmKLAAAAA0QYEFAACgCQosAAAATVBgAQAAaIICCwAAQBMUWAAAAJqgwAIAANAEBRYAAIAmKLAAAAA0QYEFAACgCQosAAAATVBgAQAAaIICCwAAQBMUWAAAAJqgwAIAANAEBRYAAIAmKLAAAAA0QYEFAACgCQosAAAATVBgAQAAaIICCwAAQBMUWAAAAJqgwAIAANAEBRYAAIAmKLAAAAA0QYEFAACgCQosAAAATVBgAQAAaIICCwAAQBMUWAAAAJqgwAIAANAEBRYAAIAmKLAAAAA0QYEFAACgCQosAAAATVBgAQAAaIICCwAAQBMUWAAAAJqgwAIAANAEBRYAAIAmKLAAAAA04dqugksp35Hk85P8VZLfTfKVtdY/3ed270gySnJvkntqrWtdzQQAAEC7unwH9pVJ/mat9ZOTvC3JN3+A2z6l1voE5RUAAICDdFZga62vqLXeM7v6m0ke2dXPAgAA4PQ7qmNgvyrJLx6wriZ5RSnljlLKTUc0DwAAAI1Z6BjYUsovJ3n4PqueW2t96ew2z01yT5L/ekDMk2ut26WUhyZ5ZSnlLbXWV+/zs25KclOSXHfddYuMDQAAQIMWKrC11s/6QOtLKV+e5POSfGattR6QsT37992llJckeVKS+xXYWuuLkrwoSdbW1vbNAgAA4PTqbBfiUsrTkjw7yRfUWv/8gNv0Sin9y5eTfHaSN3Y1EwAAAO3q8hjY70vSz3S34NeVUl6YJKWUjymlbMxu87AkF0opr0/yW0l+odb68g5nAgAAoFGdfQ5srXXlgOXbSdZnl9+e5PFdzQAAAMDpcVRnIQYAAICFKLAAAAA0QYEFAACgCQosAAAATVBgAQAAaIICCwAAQBMUWAAAAJqgwAIAANAEBRYAAIAmKLAAAAA0QYEFAACgCQosAAAATbj2uAegPaPRKMPhMJubm1ldXc1gMEi/3z/usQAAgFNOgeWqXLhwIevr65lMJhmPx+n1ernllluysbGR8+fPH/d4AADAKWYXYh6w0WiU9fX1jEajjMfjJMl4PL6yfGdn55gnBAAATjMFlgdsOBxmMpnsu24ymWQ4HB7xRAAAwFmiwPKAbW5uXnnnda/xeJytra0jnggAADhLFFgesNXV1fR6vX3X9Xq9rKysHPFEAADAWaLA8oANBoMsLe3/J7O0tJTBYHDEEwEAAGeJAssD1u/3s7GxkX6/f+Wd2F6vd2X58vLyMU8IAACcZj5Gh6ty/vz5bG9vZzgcZmtrKysrKxkMBsorAADQOQWWq7a8vJwbb7zxuMcAAADOGLsQAwAA0AQFFgAAgCYosAAAADRBgQUAAKAJCiwAAABNUGABAABoggILAABAExRYAAAAmqDAAgAA0AQFFgAAgCYosAAAADRBgQUAAKAJCiwAAABNUGABAABoggILAABAExRYAAAAmqDAAgAA0AQFFgAAgCYosAAAADRBgQUAAKAJCiwAAABNUGABAABoggILAABAExRYAAAAmqDAAgAA0AQFFgAAgCYosAAAADRBgQUAAKAJCiwAAABNUGABAABoggILAABAE6497gEAgP2NRqMMh8Nsbm5mdXU1g8Eg/X7/uMcCgGOjwALACXThwoWsr69nPB5nMpmk1+vllltuycbGRs6fP3/c4wHAsbALMQCcMKPRKOvr6xmNRplMJkmS8Xh8ZfnOzs4xTwgAx0OBBYATZjgcXimue00mkwyHwyOeCABOBgUWAE6Yzc3NjMfjfdeNx+NsbW0d8UQAcDIosABwwqyurqbX6+27rtfrZWVl5YgnAoCTobMCW0r5llLKH5RSXjf7Wj/gdk8rpby1lLJVSnlOV/MAQCsGg0GWlvbfRC8tLWUwGBzxRABwMnT9Dux31VqfMPva2LuylHJNku9P8rlJHpfkWaWUx3U8EwCcaP1+PxsbG+n3+1eKbK/Xu7J8eXn5mCcEgONx3B+j86QkW7XWtydJKeUnkzw9ye8c61QAcMzOnz+f7e3tDIfDbG1tZWVlJYPBQHkF4EzrusB+fSnly5JcTPJNtdY/2bP+XJJ37bp+V5JP63gmAGjC8vJybrzxxuMeAwBOjIV2IS6l/HIp5Y37fD09yQ8m+fgkT0hyd5L/sF/EPsvqAT/rplLKxVLKxfe85z2LjA0AAECDFnoHttb6WQ/kdqWU/5Tk5/dZdVeSR+26/sgk2wf8rBcleVGSrK2t7VtyAQAAOL26PAvxI3Zd/aIkb9znZq9JslpK+dhSyoOTPDPJy7qaCQAAgHZ1eQzs80spT8h0l+B3JPnqJCmlfEyS22ut67XWe0opX5/kl5Jck+SHa61v6nAmAAAAGtVZga21fukBy7eTrO+6vpHkfh+xAwAAALt1/TmwAAAAcCgUWAAAAJqgwAIAANAEBRYAAIAmKLAAAAA0QYEFAACgCQosAAAATVBgAQAAaIICCwAAQBMUWAAAAJqgwAIAANAEBRYAAIAmKLAAAAA0QYEFAACgCQosAAAATVBgAQAAaIICCwAAQBMUWAAAAJqgwAIAANAEBRYAAIAmKLAAAAA0QYEFAACgCQosAAAATVBgAQAAaIICCwAAQBMUWAAAAJqgwAIAANAEBRYAAIAmKLAAAAA0QYEFAACgCdce9wBAu0ajUYbDYTY3N7O6uprBYJB+v3/cYwEAcEopsMBcLly4kPX19Uwmk4zH4/R6vdxyyy3Z2NjI+fPnj3s8gFNlNBplbW0tly5dyq233uoFQ+DMKrXW457hqq2trdWLFy8e9xhwZo1Go5w7dy6j0eh+6/r9fra3t7O8vHwMkwGcPpdfMByPx5lMJun1ellaWvKCIXCqlVLuqLWu7V3uGFjgqg2Hw0wmk33XTSaTDIfDI54I4HQajUZZX1/PaDS68rg7Ho+vLN/Z2TnmCQGOlgILXLXNzc2Mx+N9143H42xtbR3xRACnkxcMAe5LgQWu2urqanq93r7rer1eVlZWjngigNPJC4YA96XAAldtMBhkaWn/h4+lpaUMBoMjngjgdPKCIcB9KbDAVev3+9nY2Ei/37/yxKrX611Z7gROAIfDC4YA9+VjdIC5nD9/Ptvb2xkOh9na2srKykoGg4HyCnCILr8weNBZiD3mAmeNj9EBADjhdnZ2vGAInCkHfYyOd2ABAE645eXl3Hjjjcc9BsCxcwwsAAAATVBgAQAAaIICCwAAQBMUWAAAAJqgwAIAANAEBRYAAIAmKLAAAAA0QYEFAACgCQosAAAATVBgAQAAaIICCwAAQBMUWAAAAJqgwAIAANAEBRYAAIAmKLAAAAA0QYEFAACgCQosAAAATVBgAQAAaMK1XQWXUoZJHju7+hFJ/rTW+oR9bveOJKMk9ya5p9a61tVMAAAAtKuzAltrHVy+XEr5D0n+7APc/Cm11j/qahYAAADa11mBvayUUpL8gyR/t+ufBQAAwOl1FMfA/q0kf1hr3TxgfU3yilLKHaWUm45gHgAAABq00DuwpZRfTvLwfVY9t9b60tnlZyX5iQ8Q8+Ra63Yp5aFJXllKeUut9dX7/KybktyUJNddd90iYwMAANCgUmvtLryUa5P8QZLra613PYDbf0uSnVrrv/9At1tbW6sXL148nCEBAAA4UUopd+x3gt+udyH+rCRvOai8llJ6pZT+5ctJPjvJGzueCQAAgAZ1XWCfmT27D5dSPqaUsjG7+rAkF0opr0/yW0l+odb68o5nAgAAoEGdnoW41voV+yzbTrI+u/z2JI/vcgYAAABOh6M4CzEAAAAsTIEFAACgCQosAAAATVBgAQAAaIICCwAAQBMUWAAAaMgNN9yQG2644bjHgGOhwAIAQCNGo1HuvvvuvP3tb8/tt9+e0Wh03CPBkVJgAQCgARcuXMi5c+eytbWVd73rXbn55ptz7ty5XLhw4bhHgyOjwAIAwAk3Go2yvr6e0WiUyWSSJBmPx1eW7+zsHPOEcDQUWAAAOOGGw+GV4rrXZDLJcDg84ongeCiwAABwwm1ubmY8Hu+7bjweZ2tr64gnguOhwAIAwAm3urqaXq+377per5eVlZUjngiOhwILAAAn3GAwyNLS/k/dl5aWMhgMjngiOB4KLAAAnHD9fj8bGxvp9/tXimyv17uyfHl5+ZgnhKNx7XEPAAAAfHDnz5/P9vZ2hsNhtra2srKyksFgoLxypiiwAADQiOXl5dx4443HPQYcG7sQAwAA0AQFFgAAgCYosAAAADRBgQUAAKAJCiwAAABNUGABAABoggILAABAExRYAAAAmqDAAgAA0AQFFgAAgCYosAAAADRBgQUAAKAJCiwAAABNUGABAABoggILAABAExRYAAAAmqDAAgAA0AQFFgAAgCYosAAAADRBgQUAAKAJCiwAAABNUGABAABoggILAABAExRYAAAAmqDAAgAA0AQFFgAAgCYosAAAAI254YYbcsMNNxz3GEfu2uMeAIC2jUajDIfDbG5uZnV1NYPBIP1+/7jHAoBTazQa5e67786lS5dy++23n6ltb6m1HvcMV21tba1evHjxuMcAOPMuXLiQ9fX1TCaTjMfj9Hq9LC0tZWNjI+fPnz/u8QDg1Lm87R2Px5lMJqd221tKuaPWuna/5QosAPMYjUY5d+5cRqPR/db1+/1sb29neXn5GCYDgNPpLG17DyqwjoEFYC7D4TCTyWTfdZPJJMPh8IgnAoDTzbZXgQVgTpubmxmPx/uuG4/H2draOuKJAOB0s+1VYAGY0+rqanq93r7rer1eVlZWjngiADjdbHsdAwvAnM7ScTgAcBKcpW2vY2ABOFT9fj8bGxvp9/tXXg3u9XpXlp+WDSgAnBS7t71LS9Mqd9a2vd6BBWAhOzs7GQ6H2draysrKSgaDwZnYgALAcdnZ2cn111+fS5cu5bbbbjuV214fowMAAEAT7EIMAABA0xRYAAAAmqDAAgAA0AQFFgAAgCYosAAAADRBgQUAAKAJCxXYUsrfL6W8qZQyKaWs7Vn3zaWUrVLKW0spn3PA939UKeWVpZTN2b8fucg8AAAAnF6LvgP7xiRfnOTVuxeWUh6X5JlJPinJ05L8QCnlmn2+/zlJfqXWuprkV2bXAQAA4H4WKrC11jfXWt+6z6qnJ/nJWutf1lp/L8lWkicdcLsXzy6/OMkXLjIPAAAAp1dXx8CeS/KuXdfvmi3b62G11ruTZPbvQzuaBwAAgMZd+8FuUEr55SQP32fVc2utLz3o2/ZZVq9msH3muCnJTbOrO6WU/d75BQAAoH2P3m/hBy2wtdbPmuOH3ZXkUbuuPzLJ9j63+8NSyiNqrXeXUh6R5N0fYI4XJXnRHLMAAABwCnS1C/HLkjyzlPIhpZSPTbKa5LcOuN2Xzy5/eZKD3tEFAADgjFv0Y3S+qJRyV5LPSPILpZRfSpJa65uS/FSS30ny8iT/pNZ67+x7bt/1kTvPS/LUUspmkqfOrgMAAMD9lFoXOjQVAAAAjkRXuxADAADAoVJgAQAAaIICCwAAQBMUWAAAAJpwKgtsKeVB+yz76I5+1icu8L2llPJppZQvnp3R+dNKKeWQ5lqbZX7+IjPuybyulPIRs8uPKaU8o5TyNxfM/OTDmG1P5oN3/x5LKU8ppXxTKeVzD/nnLJdSnnj5d3LI2V93SDlLpZSl2eUHz+b9qEPKPpL72aJ/v13ez7ivUspfP+4ZTouj/rs9jO1EF48JXT6G7foZh/J4uyfzULYPXW3Putj27so+9OcKs6wuntd8TinlB0spLyulvHR2+WmHkX3Az7t1we//nFLKjaWUx+xZ/lVz5pVSyj8opfz92eXPLKV8Tynl6y7f7w5LKeV/HELGR++5/n/N5r1pkcfH2d/VR80uP6SU8qOllDeUUoallEfOmfmdpZQnzzvTB8j9qFLKraWUfzT7f/bcUsrPl1K+o5TykQtmP6WU8n2z+8LPlFKeV0pZOYSZD/1+dqrOQlxKeUqSH0vyIUlem+SmWus7ZuvurLU+sYOf+fu11uvm+L7PTvIDSTaT/MFs8SOTrCT5ulrrK+ac5+8k+Q9J/jTJ9Ul+PclHJnlfki+ttb5rztznJPnqJH+Z5N8n+Wez7E9P8kO11u+cM/feJL+X5CeS/ESt9XfmydmT+fokN9Ra/6SU8n8n+aIkG0n+TpKLtdZvnjP3B2qtXze7fD7Jjyf53Uz/n311rXVjztxb9i5K8s1Jvi1JFvjdfmGS/5hkkuRrkvzzJOMkn5Dka2ut/33O3CO9n817H5t9b1f3s+Uk/0+SvzfL+6tM/xZeWGv9z/NkznLXknzHbNZvTvLDSZ6U5G2Z/p5fO0fmnUl+NtP71+/OO9s+uc9L8u9rrX80m/unMv1be1CSL6u1/s85c59Wa3357PKHJ/nOJJ+a5I1JvrHW+odz5j48yW2zGW9N8g2Z/v97c5J/Wmu9e87cD8v0/9Ujk/xirfXHd6278phxlZmd/N1+kJ+5yP2sk8eELh7DOny87Wr70NX27NC3vbPcQ3+u0OHzmu/O9G/pR5PcNVv8yCRflmSz1vpP58n9ID9zkfvZtyU5n+TOJJ+f5Ltrrd87WzfX/ayU8gNJHprkwUnem+l9+L8nWU/yh/P+Dkopv713Uaa/67cmSa11rhdQdv93llL+RZK/lel97fOS3FVr/cY5c3+n1vq42eVhkt9M8tNJPivJP6y1PnWOzPckeWeShyQZZno/u+pt+D65G0nekOTDkvyN2eWfyvTjSB9fa336nLnPS/KwJL+S5AszfXx4W5KvS/JttdafnjP3u9PF/azWemq+krwmySfNLj8j043/p8+uv3aB3O854Ot7k7x3zsw3J3nMPss/NsmbF5j1tUkesivrJbPLT03yigVy35TkQ5P89SSjXT+jl+SNC877N5N8a5KtJK9P8pz9fjdXkfnGXZcvJvnQ2eVrk/z2Arl37rr8q0meOLv8cZk+kZg3d5Tpg9utmT7Bvi3Jn1y+vODv9uGzv4P3JnnsbPmjF5z30O9nXdzHZrld3c9emuQrMn0QviXJv0yymuTFmT7Qz5v7W0k+N8mzkrwryTNmyz8zyW/Mmfl7mT6R/P1Z/jcm+Zh5Z9yV+4Zdl381yafOLn/Cgn9fu+9ntyf5t7O/2W9M8nML5L4809L6nCS/neTZSa6bLXvpArk/k+lnmH9hkpfNrn/I3v+WE/J329X9rKtt76E/hqW7x9uutg9dbc8Ofds7yz305wrp7nnN2w5YXjJ9Yj1v7nsP+BoluWeB3DckuXZ2+SMyfSHjuy7/jubNnP37oCR/nOTBu/6+3jBP5uz7X5bkvyT5xNn99TGZbtMeneTRi/zd7rp8Z5LervkXmfetuy7fsWfd6xaZNdPnBv9ydt94y+yx5hMWmPV19f1/p39wGLPu/lvY9f//12eXP3Le++7s+zu5n522XYgfXGt9U5LUWv9bpk8oXlxK+aIkdYHcr8z01f879nxdzPSdl3lcm/e/ErHbH2R6R5zXNbXW98wu/36mDxaptb4yybkFcu+ttV7K9BXQS5k+0KXWOl4gcxZR31hrfW6tdSXJP8701cBfK6X8v3Nmvre8f3elP0ryf8wuX5vD223+w2qtdyZJrfXtSa5ZIOuTZt/fS/IdtdZ/leRPaq3/anZ5brXW/1Vr/b0kv19rvfzK5zuz2O+hi/tZF/expLv72WNqrf+51npXnb6j8AW11s1M/zu+eIHcB9Vaf7HW+hOZ3jf+W6YXfiXv/zu+Wn9Sa/1ndfqq/zdlujG9s5Tyq6WUmxaZtZRy7ezyh9ZaXzOb9W2Zvop/GNZqrf+i1vrOWut3ZfokaF4Pq7V+b631eUk+otb672qtv1+n72A8eoHcj6+1PqfW+nO11i/I9EnV/yiL7U7d1d9tV/ezrra9XTyGdfZ4u8thbh+62p51se1Nunmu0NXzmr8opTxpn+WfmuQvFsj90ySrtdYP2/PVTzLXnh4z19Za70mSWuufZvou7IeVUn4603dQ53E5731JXlNr/avZ9XuS3DvvoLPHwp9J8qJM3xV8R5L3zR7L3zlvbpIPLaV8Sinl+kz/Lsa75p973iSvKqX861LKh84uf2FyZe+SP5szs85m26y1/pta6ycl+QeZ3ofn2iNjZmm2q/CjkiyX2e7ks23OvH8HSTIp7z8842Mye9yqtf5JpmVzXp3cz6794DdpyvtKKQ+vtf6vJKm1vqmU8plJfj7Jxy+Q+5pMX32434N6KeVb5sz84SSvKaX8ZKavSiXTP8ZnJvmhOTOT5GIp5Ycy3QXg6UleNZvzr2WxjeidpZQfz3Sj/yuZPjl5eZK/m2SRXY/uc6eotf5Wkt8qpXxTkr89Z+bXJPmvs12v3p3p7+R/JvnkzHYTm9MnznaLKUkeU0r5yDrdrWspCzyprLX+fpJnlFKenuSVpZTvWmDG+yilLNVaJ0m+ateya7LYg1wX97Mu7mNJd/ezcSnlfK31Qinl85P87ySptU5KWeg4xb+Y7T764UlqKeULa60/N9uFbpGNc2bz/VqmT1C/IdN3LwaZPsGYx/cn2ZjtdvTy2W5CP5vpu8WvW2DMh8528yyZPjkrdfZybRZ7wr77e3/0A6y7Wh+y636WWuu3llLuSvLqJMtzZnb1d9vV/ayrbe+hP4Z1+HjbyfYh3W3Putj2Jt08V+jqec1XJPnBUko/73/B6FGZvlv6FQvk/mimJXu/wx1+fJ9lD9TvllL+Tp0dnlFrvTfJjaWUf5vp4RDz+F+llOVa606t9coxibNDLhZ5USu11peUUl6R5N+UUv5RFnvecdndmR5WkiT/u5TyiFrr3bPyds8CuV+f5LmZ7eKc5BtLKeNMd6f+0jkz7/dcoNb625nuATTXrv8z357pO7nJ9HHx9lJKTfK4JIu8CPdtSV5bSnlrpu+cf20yPSY40z005vUV6eB+dtqOgf2sJO+ptb5+z/IPT/L1tdZvnTP3o5L8Ra31zw9hzN25j0vyBZm+glgy/R/7srrAsShlehKNf5zpH/Lrk/xwrfXe2atKD533la/ZOy1/P9NXlP5bpsflfUmmr4Z+/7yvrpZSvqTuOmbssMye4Hx2prszXn4345dmr1rOm7n3XZrtWuv7yvSkAn+71vqz82bv+hl/LdMHoE+rtS7yJCKllE/NdJeQv9iz/DFJztda/8ucuYd+P+vqPjbL/huZPuk5zPvZJ2e6e+snZPqO1lfVWt82e6B/Vq31e+bMfXyS52d6zN83ZroB+fJM33n7x/sVjweQ+ZO11mfOM88DyL4h0xl3389ekuRHZq+Iz5N5255FP1Brfc/sCdXza61fNmfuv559/86e5StJnldrfcacuc/PdDfGX96z/GlJvrfWujpnbhd/t11ty7ra9nbyGLYr5zAfbzvbPnS0Petq23vozxW6el6zK//h2XU/u/xCzEkz++/N7B3uvevO1Vr/4P7fNffP6mW6e+67Dynv8Uk+o9b6wsPI2yf/mkwP3Vj4sW32uHVtrfWPF8xZ3ru9OSyz/95Sa71ndp97Qqa7Ey/yDv/lbcTHJdla5PHlgOxDvZ+dqgILAGdFKeVnaq3zvvMi94gz5Xabe5iZpZRPrLW+5YPf8vTmtjRra7ktzXpYuaWUB+19YbuU8tG11j+aJ+9UHQNbSvmwUsq3l1J+rJTyJXvW/cBJyt2T+azDmvWD/MxfPOu5Lc0qt7vMLnJLKT9/mHld5rY0q9wP6OPkdpbb0qxyDz/z0M/03WBuS7O2ltvSrAvllulH89yVZLuU8opy34+Amjv3tB0D+yOZnv3wZ5J8VSnl7yX5klrrX2Z6CveTlLs38xmHMWsp5aDTqJdMdzE49bktzSq3u8wucw+wyMlEjjq3pVnlHqyrXajktjWr3DkySykHHepRMj3L71xaym1p1tZyW5q1y9xMD4v6nNm5EZ6R6bkHvrTW+puz7LmctgL78bt2H/m5UspzMz0b5BecwNyuZn1Nkv+Z/f8oPuKM5LY0q9zuMrvM3c9rDzmvy9yWZpULdOErMz0z+1/us+5Z+yw7jbktzdpabkuzdpl7n7PUl1LenORny/Qzo+d/IavO+fk7J/Er08/OW9qz7Msz/eyld56k3A5nfWOmp2/fb927zkJuS7PKbW/Wq/z5P9NKbkuzyr2S+dqOZj3zuS3NKne+zCT/I8n/ecC631tgjmZyW5q1tdyWZu0492KSh+9Z9shMP61gNG/uaXsH9r9neqr2K2eDrLW+uJTyh5l+UPtJyu1q1m/Jwcc2f8MZye0iU263uV1kdpn7QJ3V48fkHk3uszvIlNtdptxuc6828xk54HMoa60fu8AcLeW2NGtruS3N2mXuc5I8LMmVsw7XWu8q048H/Pp5Q52FGKAjpZQ7a60HHYd7onJbmvW055ZS3pAPsGtVrfWT55zhzOe2NKvc7jKv8uc3cybmrnJbmrW13JZmPUm5p+0d2Psppfx8rfXzWshtadbWcluaVW53mV3mwiG6/Pf5T2b//tjs33+YZJHPOZTb1qxyu8u8Gi3tldFVbkuztpbb0qwnJvfUF9i0dZbJlmZtLbelWeV2l9ll7n7mPsPeMeS2NOupzq21vjNJSilPrrU+edeq55RSfj3Jv55nALltzSq321mvZgS5Tc3aWm5Ls56Y3FP1ObAHaOksky3N2lpuS7PK7S6zy9z9nObjx+R2n9srpZy/fKWU8uQkvUOYRW5bs8rtLhNo0Jk8Bvak7L99XJlyu8uU223uSZnV8WNyu86dZT8x088M//DZz/izJF9Za13oRRi5bc0qt9tZH8DPfW2t9VPOcm5Ls7aW29KsJyn3LOxCvJ8Tsf/2MWbK7S5Tbre5J2VWx4/J7Sy3lHLLrqs/mum7TOPZ9adkzr0I5LY1q9zuMq/SSdor47hyW5q1tdyWZj0xuWe1wJ6I/bePMVNud5lyu809EbM6fkxux7n92b+PTfKpSV6a6btOn5/k1fPMKbfTTLnd5nYy6wPde6LW+orTmtvSrK3ltjRri7lntcACHIZeKeV8rfVCcvjHjx1ybkuznuncWuu/mmW8IskTa62j2fVvSfLT8w4ot61Z5XY7axraK6PD3JZmbS23pVnby621nrmvJK9tJbelWVvLbWlWuSdz1iRPTPL6JO9I8ntJXpfkUw5hnkPPbWlWuVcy35LkQ3Zd/5AkbzmEWc98bkuzyu181l9/IMtOc25Ls7aW29KsLeWe1XdgT8T+28eYKbe7TLnd5p6IWR0/Jrfr3JkfS/JbpZSXZLoL1hclefECeXK7zZTbbW5Xs574vTKOILelWVvLbWnWZnJP1VmIH+h+1icht6VZW8ttaVa53WV2nHvb7OLuY7JKZsdk1Vr/0UnJbWlWufvmPzHJ35pdfXU9pDOuym1rVrmdZ/5IGjgTc1e5Lc3aWm5Ls7aUe9oK7KNnF/fdz7rWOtdJOrrIbWnW1nJbmlVue7PuyX9Fkr9X339MVj/JT9dan3bScluaVS7QtT17T5Tcd++JWmv9ztOe29KsreW2NGuLuadqF+La0FkmW5q1tdyWZpXb3qx7XJfkr3Zd/6skj1kws6vclmaVC3StpTMxd5Xb0qyt5bY0a3O5p6rA7tLE/tsdZsrtLlNut7ktzZo4fkxu97lAB2pDZ2LuKrelWVvLbWnWFnNPa4H9qiQ/Ukq5z37WJzS3pVlby21pVrndZXaWW2v91lLKL+b9x2QtfJxIV7ktzSoXOEKt7ZVhD522cluatZncU1VgS0NnmWxp1tZyW5pVbneZXebuVmu9M8mdi+YcRW5Ls8oFjkhre2XYQ6et3JZmbSb3VBXYtLX/dkuztpbb0qxyu8vsMheAU6C1vTLsodNWbkuztpR7qs5CfFlp6CyTLc3aWm5Ls8ptb1YAAI7e0nEP0JEm9t/uMFNud5lyu81taVYAAI7YaduF+LIm9t/uMFNud5lyu81taVYAAI7YqdyFOElKKU/M+/ezfvVhnQ2yi9yWZm0tt6VZ5XaX2WUuAABH59QWWAAAAE6X03oMLAAAAKeMAgsAAEATFFgAAACaoMACAADQBAUWAACAJvz/rkCkvifmqKsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1152x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "constant_eventstudy = PanelOLS(df.y2, df[ind_cols_2],  \n",
    "                         entity_effects = True, time_effects = True,\n",
    "                         check_rank = False, drop_absorbed = True).fit(cov_type = 'robust')\n",
    "\n",
    "err_series = constant_eventstudy.params - constant_eventstudy.conf_int()['lower']\n",
    "cons_coef_df = pd.DataFrame({'coef': constant_eventstudy.params.values[1:],\n",
    "                        'err': err_series.values[1:],\n",
    "                        'varname': err_series.index.values[1:]\n",
    "                       })\n",
    "\n",
    "fig, ax = plt.subplots(figsize = (16, 10))\n",
    "cons_coef_df.plot(x = 'varname', y = 'coef', kind = 'bar', \n",
    "             ax = ax, color = 'none', \n",
    "             yerr ='err', legend = False)\n",
    "ax.set_ylabel('')\n",
    "ax.set_xlabel('')\n",
    "ax.set_ylim(-10, 10)\n",
    "ax.axhline(y = 0, color='r', linestyle='-')\n",
    "ax.scatter(x = np.arange(cons_coef_df.shape[0]), \n",
    "           marker = 'o', s = 50, \n",
    "           y = cons_coef_df['coef'], color = 'black')\n",
    "ax.xaxis.set_ticks_position('none')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 172,
   "id": "5c41c6ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "leads = pd.get_dummies(df['time_til'][df['time_til'] < -1], prefix = \"leads\", \n",
    "                      prefix_sep = \"_\", drop_first = False)\n",
    "\n",
    "lags = pd.get_dummies(df['time_til'][df['time_til'] > 0], prefix = \"lags\", \n",
    "                      prefix_sep = \"_\", drop_first = False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 173,
   "id": "25f13548",
   "metadata": {},
   "outputs": [],
   "source": [
    "leads_lags = pd.concat([leads, lags], axis = 0)\n",
    "df = pd.concat([df, leads_lags], axis = 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 175,
   "id": "8c11a3cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_lead_til(row):\n",
    "    if row['time_til'] < -1:\n",
    "        return row['time_til']\n",
    "    else:\n",
    "        return 0\n",
    "    \n",
    "def create_lag_til(row):\n",
    "    if row['time_til'] > 0:\n",
    "        return row['time_til']\n",
    "    else:\n",
    "        return 0\n",
    "    \n",
    "df['lead_til'] = df.apply(create_lead_til, axis = 1)\n",
    "df['lag_til'] = df.apply(create_lag_til, axis = 1)\n",
    "df['abs_lead_til'] = abs(df['lead_til'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 199,
   "id": "be106367",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "\n",
    "het_te = PanelOLS.from_formula(\"y2 ~ 1 + C(abs_lead_til):C(group) + C(lag_til):C(group) + EntityEffects\", \n",
    "                               check_rank = False, drop_absorbed = True,\n",
    "                               data = df).fit(cov_type = 'robust')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 200,
   "id": "5d95f81d",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                          PanelOLS Estimation Summary                           \n",
      "================================================================================\n",
      "Dep. Variable:                     y2   R-squared:                        0.9906\n",
      "Estimator:                   PanelOLS   R-squared (Between):             -2.4929\n",
      "No. Observations:               30000   R-squared (Within):               0.9906\n",
      "Date:                Tue, Feb 15 2022   R-squared (Overall):              0.7486\n",
      "Time:                        22:11:56   Log-likelihood                -4.538e+04\n",
      "Cov. Estimator:                Robust                                           \n",
      "                                        F-statistic:                   2.704e+04\n",
      "Entities:                        1000   P-value                           0.0000\n",
      "Avg Obs:                       30.000   Distribution:               F(112,28888)\n",
      "Min Obs:                       30.000                                           \n",
      "Max Obs:                       30.000   F-statistic (robust):          5.191e+05\n",
      "                                        P-value                           0.0000\n",
      "Time periods:                      30   Distribution:               F(112,28888)\n",
      "Avg Obs:                      1000.00                                           \n",
      "Min Obs:                      1000.00                                           \n",
      "Max Obs:                      1000.00                                           \n",
      "                                                                                \n",
      "                                           Parameter Estimates                                           \n",
      "=========================================================================================================\n",
      "                                       Parameter  Std. Err.     T-stat    P-value    Lower CI    Upper CI\n",
      "---------------------------------------------------------------------------------------------------------\n",
      "Intercept                                 21.408     0.0946     226.30     0.0000      21.222      21.593\n",
      "C(abs_lead_til)[T.2.0]:C(group)[T.1]     -6.5210     0.2516    -25.922     0.0000     -7.0141     -6.0279\n",
      "C(abs_lead_til)[T.2.0]:C(group)[T.2]     -5.4918     0.2051    -26.779     0.0000     -5.8938     -5.0899\n",
      "C(abs_lead_til)[T.2.0]:C(group)[T.3]     -4.5003     0.1597    -28.174     0.0000     -4.8134     -4.1872\n",
      "C(abs_lead_til)[T.2.0]:C(group)[T.4]     -3.5116     0.1154    -30.441     0.0000     -3.7377     -3.2855\n",
      "C(abs_lead_til)[T.3.0]:C(group)[T.1]     -7.5084     0.2515    -29.849     0.0000     -8.0015     -7.0154\n",
      "C(abs_lead_til)[T.3.0]:C(group)[T.2]     -6.4766     0.2052    -31.569     0.0000     -6.8787     -6.0745\n",
      "C(abs_lead_til)[T.3.0]:C(group)[T.3]     -5.4964     0.1597    -34.410     0.0000     -5.8095     -5.1833\n",
      "C(abs_lead_til)[T.3.0]:C(group)[T.4]     -4.4931     0.1153    -38.966     0.0000     -4.7191     -4.2671\n",
      "C(abs_lead_til)[T.4.0]:C(group)[T.1]     -8.4830     0.2515    -33.724     0.0000     -8.9760     -7.9900\n",
      "C(abs_lead_til)[T.4.0]:C(group)[T.2]     -7.4989     0.2051    -36.558     0.0000     -7.9009     -7.0968\n",
      "C(abs_lead_til)[T.4.0]:C(group)[T.3]     -6.4816     0.1599    -40.533     0.0000     -6.7950     -6.1682\n",
      "C(abs_lead_til)[T.4.0]:C(group)[T.4]     -5.5037     0.1153    -47.726     0.0000     -5.7297     -5.2776\n",
      "C(abs_lead_til)[T.5.0]:C(group)[T.1]     -9.4817     0.2516    -37.687     0.0000     -9.9748     -8.9886\n",
      "C(abs_lead_til)[T.5.0]:C(group)[T.2]     -8.4806     0.2050    -41.360     0.0000     -8.8825     -8.0787\n",
      "C(abs_lead_til)[T.5.0]:C(group)[T.3]     -7.4919     0.1598    -46.897     0.0000     -7.8051     -7.1788\n",
      "C(abs_lead_til)[T.5.0]:C(group)[T.4]     -6.4925     0.1153    -56.328     0.0000     -6.7184     -6.2666\n",
      "C(abs_lead_til)[T.6.0]:C(group)[T.1]     -10.494     0.2516    -41.710     0.0000     -10.987     -10.001\n",
      "C(abs_lead_til)[T.6.0]:C(group)[T.2]     -9.4829     0.2052    -46.207     0.0000     -9.8852     -9.0807\n",
      "C(abs_lead_til)[T.6.0]:C(group)[T.3]     -8.4779     0.1599    -53.034     0.0000     -8.7912     -8.1645\n",
      "C(abs_lead_til)[T.6.0]:C(group)[T.4]     -7.4927     0.1153    -64.983     0.0000     -7.7187     -7.2667\n",
      "C(abs_lead_til)[T.7.0]:C(group)[T.2]     -10.500     0.2051    -51.187     0.0000     -10.902     -10.098\n",
      "C(abs_lead_til)[T.7.0]:C(group)[T.3]     -9.4903     0.1597    -59.424     0.0000     -9.8034     -9.1773\n",
      "C(abs_lead_til)[T.7.0]:C(group)[T.4]     -8.5021     0.1151    -73.846     0.0000     -8.7277     -8.2764\n",
      "C(abs_lead_til)[T.8.0]:C(group)[T.2]     -11.479     0.2051    -55.972     0.0000     -11.881     -11.077\n",
      "C(abs_lead_til)[T.8.0]:C(group)[T.3]     -10.529     0.1597    -65.931     0.0000     -10.842     -10.216\n",
      "C(abs_lead_til)[T.8.0]:C(group)[T.4]     -9.5159     0.1152    -82.569     0.0000     -9.7418     -9.2900\n",
      "C(abs_lead_til)[T.9.0]:C(group)[T.2]     -12.488     0.2052    -60.859     0.0000     -12.890     -12.085\n",
      "C(abs_lead_til)[T.9.0]:C(group)[T.3]     -11.529     0.1598    -72.150     0.0000     -11.842     -11.215\n",
      "C(abs_lead_til)[T.9.0]:C(group)[T.4]     -10.501     0.1154    -91.004     0.0000     -10.727     -10.275\n",
      "C(abs_lead_til)[T.10.0]:C(group)[T.2]    -13.469     0.2051    -65.685     0.0000     -13.871     -13.067\n",
      "C(abs_lead_til)[T.10.0]:C(group)[T.3]    -12.498     0.1599    -78.184     0.0000     -12.811     -12.184\n",
      "C(abs_lead_til)[T.10.0]:C(group)[T.4]    -11.508     0.1154    -99.769     0.0000     -11.735     -11.282\n",
      "C(abs_lead_til)[T.11.0]:C(group)[T.2]    -14.478     0.2052    -70.570     0.0000     -14.880     -14.076\n",
      "C(abs_lead_til)[T.11.0]:C(group)[T.3]    -13.521     0.1597    -84.687     0.0000     -13.834     -13.208\n",
      "C(abs_lead_til)[T.11.0]:C(group)[T.4]    -12.495     0.1153    -108.37     0.0000     -12.721     -12.269\n",
      "C(abs_lead_til)[T.12.0]:C(group)[T.2]    -15.484     0.2052    -75.471     0.0000     -15.886     -15.082\n",
      "C(abs_lead_til)[T.12.0]:C(group)[T.3]    -14.503     0.1598    -90.773     0.0000     -14.816     -14.190\n",
      "C(abs_lead_til)[T.12.0]:C(group)[T.4]    -13.494     0.1153    -117.05     0.0000     -13.720     -13.268\n",
      "C(abs_lead_til)[T.13.0]:C(group)[T.3]    -15.508     0.1597    -97.135     0.0000     -15.821     -15.195\n",
      "C(abs_lead_til)[T.13.0]:C(group)[T.4]    -14.505     0.1154    -125.71     0.0000     -14.731     -14.279\n",
      "C(abs_lead_til)[T.14.0]:C(group)[T.3]    -16.496     0.1598    -103.21     0.0000     -16.810     -16.183\n",
      "C(abs_lead_til)[T.14.0]:C(group)[T.4]    -15.498     0.1154    -134.34     0.0000     -15.724     -15.272\n",
      "C(abs_lead_til)[T.15.0]:C(group)[T.3]    -17.501     0.1597    -109.56     0.0000     -17.814     -17.188\n",
      "C(abs_lead_til)[T.15.0]:C(group)[T.4]    -16.502     0.1153    -143.17     0.0000     -16.728     -16.276\n",
      "C(abs_lead_til)[T.16.0]:C(group)[T.3]    -18.488     0.1597    -115.78     0.0000     -18.801     -18.175\n",
      "C(abs_lead_til)[T.16.0]:C(group)[T.4]    -17.513     0.1153    -151.95     0.0000     -17.739     -17.287\n",
      "C(abs_lead_til)[T.17.0]:C(group)[T.3]    -19.515     0.1598    -122.15     0.0000     -19.828     -19.201\n",
      "C(abs_lead_til)[T.17.0]:C(group)[T.4]    -18.481     0.1153    -160.32     0.0000     -18.707     -18.255\n",
      "C(abs_lead_til)[T.18.0]:C(group)[T.3]    -20.483     0.1598    -128.22     0.0000     -20.796     -20.170\n",
      "C(abs_lead_til)[T.18.0]:C(group)[T.4]    -19.525     0.1153    -169.30     0.0000     -19.751     -19.299\n",
      "C(abs_lead_til)[T.19.0]:C(group)[T.4]    -20.506     0.1154    -177.66     0.0000     -20.732     -20.280\n",
      "C(abs_lead_til)[T.20.0]:C(group)[T.4]    -21.527     0.1154    -186.53     0.0000     -21.753     -21.300\n",
      "C(abs_lead_til)[T.21.0]:C(group)[T.4]    -22.497     0.1153    -195.12     0.0000     -22.723     -22.271\n",
      "C(abs_lead_til)[T.22.0]:C(group)[T.4]    -23.483     0.1151    -203.98     0.0000     -23.709     -23.258\n",
      "C(abs_lead_til)[T.23.0]:C(group)[T.4]    -24.476     0.1153    -212.31     0.0000     -24.702     -24.250\n",
      "C(abs_lead_til)[T.24.0]:C(group)[T.4]    -25.510     0.1154    -221.13     0.0000     -25.736     -25.284\n",
      "C(group)[T.1]:C(lag_til)[T.1.0]           6.4864     0.2516     25.776     0.0000      5.9932      6.9796\n",
      "C(group)[T.1]:C(lag_til)[T.2.0]           7.4991     0.2516     29.807     0.0000      7.0059      7.9922\n",
      "C(group)[T.1]:C(lag_til)[T.3.0]           8.4780     0.2516     33.694     0.0000      7.9848      8.9712\n",
      "C(group)[T.1]:C(lag_til)[T.4.0]           9.5190     0.2516     37.832     0.0000      9.0258      10.012\n",
      "C(group)[T.1]:C(lag_til)[T.5.0]           10.516     0.2515     41.813     0.0000      10.023      11.009\n",
      "C(group)[T.1]:C(lag_til)[T.6.0]           11.515     0.2516     45.758     0.0000      11.021      12.008\n",
      "C(group)[T.1]:C(lag_til)[T.7.0]           12.486     0.2516     49.626     0.0000      11.992      12.979\n",
      "C(group)[T.1]:C(lag_til)[T.8.0]           13.519     0.2516     53.725     0.0000      13.026      14.012\n",
      "C(group)[T.1]:C(lag_til)[T.9.0]           14.496     0.2515     57.628     0.0000      14.003      14.989\n",
      "C(group)[T.1]:C(lag_til)[T.10.0]          15.527     0.2515     61.726     0.0000      15.034      16.020\n",
      "C(group)[T.1]:C(lag_til)[T.11.0]          16.505     0.2515     65.622     0.0000      16.012      16.997\n",
      "C(group)[T.1]:C(lag_til)[T.12.0]          17.487     0.2516     69.506     0.0000      16.994      17.980\n",
      "C(group)[T.1]:C(lag_til)[T.13.0]          18.503     0.2516     73.545     0.0000      18.010      18.996\n",
      "C(group)[T.1]:C(lag_til)[T.14.0]          19.482     0.2516     77.429     0.0000      18.989      19.975\n",
      "C(group)[T.1]:C(lag_til)[T.15.0]          20.485     0.2516     81.419     0.0000      19.991      20.978\n",
      "C(group)[T.1]:C(lag_til)[T.16.0]          21.492     0.2516     85.437     0.0000      20.999      21.985\n",
      "C(group)[T.1]:C(lag_til)[T.17.0]          22.519     0.2516     89.502     0.0000      22.026      23.012\n",
      "C(group)[T.1]:C(lag_til)[T.18.0]          23.501     0.2516     93.407     0.0000      23.008      23.995\n",
      "C(group)[T.1]:C(lag_til)[T.19.0]          24.508     0.2516     97.422     0.0000      24.015      25.001\n",
      "C(group)[T.1]:C(lag_til)[T.20.0]          25.508     0.2515     101.41     0.0000      25.015      26.001\n",
      "C(group)[T.1]:C(lag_til)[T.21.0]          26.516     0.2516     105.38     0.0000      26.022      27.009\n",
      "C(group)[T.1]:C(lag_til)[T.22.0]          27.517     0.2516     109.36     0.0000      27.024      28.010\n",
      "C(group)[T.1]:C(lag_til)[T.23.0]          28.500     0.2516     113.29     0.0000      28.007      28.993\n",
      "C(group)[T.2]:C(lag_til)[T.1.0]           5.5187     0.2052     26.896     0.0000      5.1165      5.9209\n",
      "C(group)[T.2]:C(lag_til)[T.2.0]           6.5243     0.2051     31.809     0.0000      6.1223      6.9263\n",
      "C(group)[T.2]:C(lag_til)[T.3.0]           7.5184     0.2052     36.637     0.0000      7.1162      7.9206\n",
      "C(group)[T.2]:C(lag_til)[T.4.0]           8.5301     0.2051     41.589     0.0000      8.1281      8.9321\n",
      "C(group)[T.2]:C(lag_til)[T.5.0]           9.5271     0.2051     46.440     0.0000      9.1250      9.9293\n",
      "C(group)[T.2]:C(lag_til)[T.6.0]           10.512     0.2052     51.235     0.0000      10.110      10.914\n",
      "C(group)[T.2]:C(lag_til)[T.7.0]           11.515     0.2051     56.154     0.0000      11.113      11.917\n",
      "C(group)[T.2]:C(lag_til)[T.8.0]           12.537     0.2052     61.104     0.0000      12.134      12.939\n",
      "C(group)[T.2]:C(lag_til)[T.9.0]           13.508     0.2051     65.869     0.0000      13.106      13.910\n",
      "C(group)[T.2]:C(lag_til)[T.10.0]          14.517     0.2052     70.749     0.0000      14.115      14.919\n",
      "C(group)[T.2]:C(lag_til)[T.11.0]          15.520     0.2052     75.627     0.0000      15.118      15.922\n",
      "C(group)[T.2]:C(lag_til)[T.12.0]          16.506     0.2051     80.464     0.0000      16.104      16.908\n",
      "C(group)[T.2]:C(lag_til)[T.13.0]          17.507     0.2051     85.359     0.0000      17.105      17.909\n",
      "C(group)[T.2]:C(lag_til)[T.14.0]          18.538     0.2051     90.373     0.0000      18.135      18.940\n",
      "C(group)[T.2]:C(lag_til)[T.15.0]          19.498     0.2052     95.017     0.0000      19.096      19.900\n",
      "C(group)[T.2]:C(lag_til)[T.16.0]          20.539     0.2052     100.09     0.0000      20.137      20.942\n",
      "C(group)[T.2]:C(lag_til)[T.17.0]          21.524     0.2052     104.88     0.0000      21.121      21.926\n",
      "C(group)[T.3]:C(lag_til)[T.1.0]           4.5031     0.1598     28.188     0.0000      4.1900      4.8163\n",
      "C(group)[T.3]:C(lag_til)[T.2.0]           5.4955     0.1597     34.405     0.0000      5.1824      5.8086\n",
      "C(group)[T.3]:C(lag_til)[T.3.0]           6.5162     0.1597     40.799     0.0000      6.2031      6.8292\n",
      "C(group)[T.3]:C(lag_til)[T.4.0]           7.4803     0.1598     46.816     0.0000      7.1672      7.7935\n",
      "C(group)[T.3]:C(lag_til)[T.5.0]           8.4903     0.1598     53.144     0.0000      8.1772      8.8035\n",
      "C(group)[T.3]:C(lag_til)[T.6.0]           9.5041     0.1598     59.481     0.0000      9.1909      9.8173\n",
      "C(group)[T.3]:C(lag_til)[T.7.0]           10.511     0.1597     65.819     0.0000      10.198      10.824\n",
      "C(group)[T.3]:C(lag_til)[T.8.0]           11.498     0.1597     72.006     0.0000      11.185      11.811\n",
      "C(group)[T.3]:C(lag_til)[T.9.0]           12.493     0.1597     78.229     0.0000      12.180      12.806\n",
      "C(group)[T.3]:C(lag_til)[T.10.0]          13.511     0.1598     84.555     0.0000      13.197      13.824\n",
      "C(group)[T.3]:C(lag_til)[T.11.0]          14.489     0.1597     90.753     0.0000      14.176      14.802\n",
      "C(group)[T.4]:C(lag_til)[T.1.0]           3.4826     0.1155     30.157     0.0000      3.2562      3.7089\n",
      "C(group)[T.4]:C(lag_til)[T.2.0]           4.5105     0.1153     39.133     0.0000      4.2846      4.7364\n",
      "C(group)[T.4]:C(lag_til)[T.3.0]           5.4688     0.1152     47.470     0.0000      5.2430      5.6946\n",
      "C(group)[T.4]:C(lag_til)[T.4.0]           6.4992     0.1154     56.331     0.0000      6.2730      6.7253\n",
      "C(group)[T.4]:C(lag_til)[T.5.0]           7.5052     0.1153     65.093     0.0000      7.2792      7.7312\n",
      "=========================================================================================================\n",
      "\n",
      "F-test for Poolability: 100.19\n",
      "P-value: 0.0000\n",
      "Distribution: F(999,28888)\n",
      "\n",
      "Included effects: Entity\n"
     ]
    }
   ],
   "source": [
    "print(het_te)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 177,
   "id": "e09431c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "leads_3 = {'dd_-2', 'dd_-1'}\n",
    "ind_cols_3 = [ind for ind in dd if ind not in leads_3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 202,
   "id": "cb4fb958",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                          PanelOLS Estimation Summary                           \n",
      "================================================================================\n",
      "Dep. Variable:                     y2   R-squared:                        0.9380\n",
      "Estimator:                   PanelOLS   R-squared (Between):             -0.7056\n",
      "No. Observations:               30000   R-squared (Within):              -0.1219\n",
      "Date:                Tue, Feb 15 2022   R-squared (Overall):             -0.1624\n",
      "Time:                        22:13:21   Log-likelihood                -2.308e+04\n",
      "Cov. Estimator:                Robust                                           \n",
      "                                        F-statistic:                   1.067e+04\n",
      "Entities:                        1000   P-value                           0.0000\n",
      "Avg Obs:                       30.000   Distribution:                F(41,28930)\n",
      "Min Obs:                       30.000                                           \n",
      "Max Obs:                       30.000   F-statistic (robust):             8822.8\n",
      "                                        P-value                           0.0000\n",
      "Time periods:                      30   Distribution:                F(41,28930)\n",
      "Avg Obs:                      1000.00                                           \n",
      "Min Obs:                      1000.00                                           \n",
      "Max Obs:                      1000.00                                           \n",
      "                                                                                \n",
      "                             Parameter Estimates                              \n",
      "==============================================================================\n",
      "            Parameter  Std. Err.     T-stat    P-value    Lower CI    Upper CI\n",
      "------------------------------------------------------------------------------\n",
      "dd_-24         9.4353     0.0909     103.83     0.0000      9.2572      9.6135\n",
      "dd_-23         9.4691     0.0905     104.62     0.0000      9.2917      9.6465\n",
      "dd_-22         9.4477     0.0902     104.75     0.0000      9.2709      9.6245\n",
      "dd_-21         9.4354     0.0907     103.98     0.0000      9.2575      9.6132\n",
      "dd_-20         9.4296     0.0478     197.27     0.0000      9.3359      9.5233\n",
      "dd_-19         9.4480     0.0513     184.20     0.0000      9.3474      9.5485\n",
      "dd_-18         6.1454     0.0818     75.085     0.0000      5.9849      6.3058\n",
      "dd_-17         6.1484     0.0813     75.647     0.0000      5.9891      6.3078\n",
      "dd_-16         6.1412     0.0813     75.564     0.0000      5.9819      6.3005\n",
      "dd_-15         6.1328     0.0816     75.189     0.0000      5.9730      6.2927\n",
      "dd_-14         6.1600     0.0396     155.62     0.0000      6.0824      6.2376\n",
      "dd_-13         6.1398     0.0422     145.46     0.0000      6.0571      6.2226\n",
      "dd_-12         2.9998     0.0674     44.513     0.0000      2.8677      3.1319\n",
      "dd_-11         2.9904     0.0673     44.433     0.0000      2.8585      3.1223\n",
      "dd_-10         2.9939     0.0671     44.634     0.0000      2.8624      3.1254\n",
      "dd_-9          2.9686     0.0675     43.981     0.0000      2.8363      3.1009\n",
      "dd_-8          2.9914     0.0256     116.70     0.0000      2.9412      3.0417\n",
      "dd_-7          2.9899     0.0289     103.51     0.0000      2.9333      3.0465\n",
      "dd_-6          0.0098     0.0610     0.1612     0.8719     -0.1098      0.1294\n",
      "dd_-5          0.0016     0.0608     0.0261     0.9792     -0.1176      0.1208\n",
      "dd_-4         -0.0020     0.0607    -0.0322     0.9743     -0.1209      0.1170\n",
      "dd_-3         -0.0206     0.0610    -0.3381     0.7353     -0.1402      0.0990\n",
      "dd_0           4.6371     0.0454     102.06     0.0000      4.5481      4.7262\n",
      "dd_1           4.6354     0.0457     101.32     0.0000      4.5457      4.7250\n",
      "dd_2           4.6501     0.0452     102.94     0.0000      4.5616      4.7387\n",
      "dd_3           4.6190     0.0459     100.73     0.0000      4.5291      4.7089\n",
      "dd_4           4.6648     0.0383     121.86     0.0000      4.5898      4.7398\n",
      "dd_5           4.6484     0.0301     154.55     0.0000      4.5895      4.7074\n",
      "dd_6           2.9556     0.0330     89.440     0.0000      2.8908      3.0204\n",
      "dd_7           2.9380     0.0329     89.284     0.0000      2.8735      3.0025\n",
      "dd_8           2.9571     0.0323     91.502     0.0000      2.8938      3.0204\n",
      "dd_9           2.9181     0.0332     87.903     0.0000      2.8530      2.9831\n",
      "dd_10          2.9755     0.0402     74.058     0.0000      2.8967      3.0542\n",
      "dd_11          2.9429     0.0201     146.61     0.0000      2.9035      2.9822\n",
      "dd_12          1.3897     0.0318     43.736     0.0000      1.3274      1.4520\n",
      "dd_13          1.3856     0.0314     44.102     0.0000      1.3240      1.4472\n",
      "dd_14          1.3982     0.0317     44.100     0.0000      1.3361      1.4604\n",
      "dd_15          1.3594     0.0323     42.046     0.0000      1.2960      1.4228\n",
      "dd_16          1.4259     0.0553     25.788     0.0000      1.3175      1.5342\n",
      "dd_17          1.4106     0.0285     49.443     0.0000      1.3547      1.4665\n",
      "dd_22          0.0284     0.0613     0.4634     0.6431     -0.0917      0.1485\n",
      "cons           18.867     0.0285     662.38     0.0000      18.811      18.922\n",
      "==============================================================================\n",
      "\n",
      "F-test for Poolability: 3583.0\n",
      "P-value: 0.0000\n",
      "Distribution: F(1028,28930)\n",
      "\n",
      "Included effects: Entity, Time\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\sajidmu2\\Anaconda3\\lib\\site-packages\\linearmodels\\panel\\model.py:1831: AbsorbingEffectWarning: \n",
      "Variables have been fully absorbed and have removed from the regression:\n",
      "\n",
      "dd_18, dd_19, dd_20, dd_21, dd_23\n",
      "\n",
      "  warnings.warn(\n"
     ]
    }
   ],
   "source": [
    "het_te = PanelOLS(df.y2, df[ind_cols_3 + ['cons']], entity_effects = True, time_effects = True, \n",
    "                  check_rank = False, drop_absorbed = True).fit(cov_type = 'robust')\n",
    "print(het_te)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "caac1d51",
   "metadata": {},
   "source": [
    "### Bacon decomposition shows the problem -- notice all those late to early 2x2s!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 203,
   "id": "c5c4b5ef",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<rpy2.rinterface_lib.sexp.NULLType object at 0x00000200D4FD3500> [RTYPES.NILSXP]"
      ]
     },
     "execution_count": 203,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import rpy2.robjects as ro\n",
    "from rpy2.robjects.packages import importr\n",
    "from rpy2.robjects import pandas2ri\n",
    "from rpy2.robjects.conversion import localconverter\n",
    "from rpy2.robjects import IntVector, Formula\n",
    "import rpy2.robjects.packages as rpackages\n",
    "utils = rpackages.importr('utils')\n",
    "utils.chooseCRANmirror(ind = 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 204,
   "id": "3ed58fce",
   "metadata": {},
   "outputs": [],
   "source": [
    "with localconverter(ro.default_converter + pandas2ri.converter):\n",
    "      rdf = ro.conversion.py2rpy(df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 205,
   "id": "92e948ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "\n",
    "# Bacon Decomposition\n",
    "utils.install_packages(\"bacondecomp\")\n",
    "bacondecomp = rpackages.importr('bacondecomp')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 209,
   "id": "ff8450c3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                      type weight   avg_est\n",
      "1 Earlier vs Later Treated    0.5  51.79904\n",
      "2 Later vs Earlier Treated    0.5 -65.20756\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>treated</th>\n",
       "      <th>untreated</th>\n",
       "      <th>estimate</th>\n",
       "      <th>weight</th>\n",
       "      <th>type</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1992</td>\n",
       "      <td>1986</td>\n",
       "      <td>-44.002869</td>\n",
       "      <td>0.100000</td>\n",
       "      <td>Later vs Earlier Treated</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1998</td>\n",
       "      <td>1986</td>\n",
       "      <td>-81.012049</td>\n",
       "      <td>0.133333</td>\n",
       "      <td>Later vs Earlier Treated</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2004</td>\n",
       "      <td>1986</td>\n",
       "      <td>-106.007429</td>\n",
       "      <td>0.100000</td>\n",
       "      <td>Later vs Earlier Treated</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>1986</td>\n",
       "      <td>1992</td>\n",
       "      <td>34.999866</td>\n",
       "      <td>0.033333</td>\n",
       "      <td>Earlier vs Later Treated</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>1998</td>\n",
       "      <td>1992</td>\n",
       "      <td>-33.015604</td>\n",
       "      <td>0.066667</td>\n",
       "      <td>Later vs Earlier Treated</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>2004</td>\n",
       "      <td>1992</td>\n",
       "      <td>-57.995785</td>\n",
       "      <td>0.066667</td>\n",
       "      <td>Later vs Earlier Treated</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>1986</td>\n",
       "      <td>1998</td>\n",
       "      <td>65.000634</td>\n",
       "      <td>0.066667</td>\n",
       "      <td>Earlier vs Later Treated</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>1992</td>\n",
       "      <td>1998</td>\n",
       "      <td>27.987025</td>\n",
       "      <td>0.066667</td>\n",
       "      <td>Earlier vs Later Treated</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>2004</td>\n",
       "      <td>1998</td>\n",
       "      <td>-22.011540</td>\n",
       "      <td>0.033333</td>\n",
       "      <td>Later vs Earlier Treated</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>1986</td>\n",
       "      <td>2004</td>\n",
       "      <td>95.004862</td>\n",
       "      <td>0.100000</td>\n",
       "      <td>Earlier vs Later Treated</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>1992</td>\n",
       "      <td>2004</td>\n",
       "      <td>52.006278</td>\n",
       "      <td>0.133333</td>\n",
       "      <td>Earlier vs Later Treated</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>1998</td>\n",
       "      <td>2004</td>\n",
       "      <td>20.990225</td>\n",
       "      <td>0.100000</td>\n",
       "      <td>Earlier vs Later Treated</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    treated  untreated    estimate    weight                      type\n",
       "2      1992       1986  -44.002869  0.100000  Later vs Earlier Treated\n",
       "3      1998       1986  -81.012049  0.133333  Later vs Earlier Treated\n",
       "4      2004       1986 -106.007429  0.100000  Later vs Earlier Treated\n",
       "5      1986       1992   34.999866  0.033333  Earlier vs Later Treated\n",
       "7      1998       1992  -33.015604  0.066667  Later vs Earlier Treated\n",
       "8      2004       1992  -57.995785  0.066667  Later vs Earlier Treated\n",
       "9      1986       1998   65.000634  0.066667  Earlier vs Later Treated\n",
       "10     1992       1998   27.987025  0.066667  Earlier vs Later Treated\n",
       "12     2004       1998  -22.011540  0.033333  Later vs Earlier Treated\n",
       "13     1986       2004   95.004862  0.100000  Earlier vs Later Treated\n",
       "14     1992       2004   52.006278  0.133333  Earlier vs Later Treated\n",
       "15     1998       2004   20.990225  0.100000  Earlier vs Later Treated"
      ]
     },
     "execution_count": 209,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "get_df_y = bacondecomp.bacon(Formula('y ~ treat'), data = rdf, \n",
    "                  id_var = 'id', \n",
    "                  time_var = 'year')\n",
    "\n",
    "with localconverter(ro.default_converter + pandas2ri.converter):\n",
    "      decomp_df_y = ro.conversion.rpy2py(get_df_y)\n",
    "        \n",
    "decomp_df_y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 210,
   "id": "201a41ae",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                      type weight avg_est\n",
      "1 Earlier vs Later Treated    0.5 8.39645\n",
      "2 Later vs Earlier Treated    0.5 5.59325\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>treated</th>\n",
       "      <th>untreated</th>\n",
       "      <th>estimate</th>\n",
       "      <th>weight</th>\n",
       "      <th>type</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1992</td>\n",
       "      <td>1986</td>\n",
       "      <td>8.000033</td>\n",
       "      <td>0.100000</td>\n",
       "      <td>Later vs Earlier Treated</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1998</td>\n",
       "      <td>1986</td>\n",
       "      <td>5.999868</td>\n",
       "      <td>0.133333</td>\n",
       "      <td>Later vs Earlier Treated</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2004</td>\n",
       "      <td>1986</td>\n",
       "      <td>3.987514</td>\n",
       "      <td>0.100000</td>\n",
       "      <td>Later vs Earlier Treated</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>1986</td>\n",
       "      <td>1992</td>\n",
       "      <td>10.001493</td>\n",
       "      <td>0.033333</td>\n",
       "      <td>Earlier vs Later Treated</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>1998</td>\n",
       "      <td>1992</td>\n",
       "      <td>5.983675</td>\n",
       "      <td>0.066667</td>\n",
       "      <td>Later vs Earlier Treated</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>2004</td>\n",
       "      <td>1992</td>\n",
       "      <td>3.992517</td>\n",
       "      <td>0.066667</td>\n",
       "      <td>Later vs Earlier Treated</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>1986</td>\n",
       "      <td>1998</td>\n",
       "      <td>10.003778</td>\n",
       "      <td>0.066667</td>\n",
       "      <td>Earlier vs Later Treated</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>1992</td>\n",
       "      <td>1998</td>\n",
       "      <td>7.983525</td>\n",
       "      <td>0.066667</td>\n",
       "      <td>Earlier vs Later Treated</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>2004</td>\n",
       "      <td>1998</td>\n",
       "      <td>3.984185</td>\n",
       "      <td>0.033333</td>\n",
       "      <td>Later vs Earlier Treated</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>1986</td>\n",
       "      <td>2004</td>\n",
       "      <td>10.001248</td>\n",
       "      <td>0.100000</td>\n",
       "      <td>Earlier vs Later Treated</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>1992</td>\n",
       "      <td>2004</td>\n",
       "      <td>7.999977</td>\n",
       "      <td>0.133333</td>\n",
       "      <td>Earlier vs Later Treated</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>1998</td>\n",
       "      <td>2004</td>\n",
       "      <td>5.989015</td>\n",
       "      <td>0.100000</td>\n",
       "      <td>Earlier vs Later Treated</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    treated  untreated   estimate    weight                      type\n",
       "2      1992       1986   8.000033  0.100000  Later vs Earlier Treated\n",
       "3      1998       1986   5.999868  0.133333  Later vs Earlier Treated\n",
       "4      2004       1986   3.987514  0.100000  Later vs Earlier Treated\n",
       "5      1986       1992  10.001493  0.033333  Earlier vs Later Treated\n",
       "7      1998       1992   5.983675  0.066667  Later vs Earlier Treated\n",
       "8      2004       1992   3.992517  0.066667  Later vs Earlier Treated\n",
       "9      1986       1998  10.003778  0.066667  Earlier vs Later Treated\n",
       "10     1992       1998   7.983525  0.066667  Earlier vs Later Treated\n",
       "12     2004       1998   3.984185  0.033333  Later vs Earlier Treated\n",
       "13     1986       2004  10.001248  0.100000  Earlier vs Later Treated\n",
       "14     1992       2004   7.999977  0.133333  Earlier vs Later Treated\n",
       "15     1998       2004   5.989015  0.100000  Earlier vs Later Treated"
      ]
     },
     "execution_count": 210,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "get_df_y2 = bacondecomp.bacon(Formula('y2 ~ treat'), data = rdf, \n",
    "                  id_var = 'id', \n",
    "                  time_var = 'year')\n",
    "\n",
    "with localconverter(ro.default_converter + pandas2ri.converter):\n",
    "      decomp_df_y2 = ro.conversion.rpy2py(get_df_y2)\n",
    "        \n",
    "decomp_df_y2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d0e1fceb",
   "metadata": {},
   "source": [
    "* Callaway and Sant'anna (2018)\n",
    "* When there are no coariates X and everybody eventually gets treated, Callaway & Sant'Anna (2018) suggest to estimate ATT(g,t) using\n",
    "* ATT(g,t) = E [ (Gg/E[Gg] - ( 1-Dt)p/1-p/E[1-Dt] )( Yt - Y_{g-1}) ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 181,
   "id": "a45c0c38",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create group indicators\n",
    "df = pd.concat([df, pd.get_dummies(df['group'], prefix = 'g', prefix_sep = \"\")], axis = 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad7aefbc",
   "metadata": {},
   "source": [
    "* Estimate propensity scores. \n",
    "* In this setup where there is no \"never treated\" group, we need propensity scores per time period. \n",
    "* We can exploit the staggered rollout to simplify this a bit!\n",
    "* Since I only know how to work things out well in the \"wide\" data format instead of the \"long\" data format, I simplify things a bit!\n",
    "* pscore for Group 1986 that can be used until year 1991 (since at 1992 another set of units get treated)\n",
    "* Note that I am using only a single year of data (this is because I want to use \"wide format\".)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 182,
   "id": "48d9c8e0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.562335\n",
      "         Iterations 5\n"
     ]
    }
   ],
   "source": [
    "logit_mod = sm.Logit.from_formula('g1 ~ 1', \n",
    "                                  data = df.query('(year == 1991) & ((g1 == 1) | (time_til < 0))')).fit()\n",
    "df['pg1_1991'] = logit_mod.predict(df['g1']).g1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44803396",
   "metadata": {},
   "source": [
    "* Calculate ATT(1986,1986)\n",
    "* This the the ATT for the group first treated at 1986 (g1==1) in the first period since treatment (1986)\n",
    "* Gen outcomes (easiest way to transform data into wide form is this)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 217,
   "id": "64ab77c3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10.006858432324133"
      ]
     },
     "execution_count": 217,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df['ypost'] = df['y'][df['year'] == 1986] \n",
    "\n",
    "df['ypre'] = df['y'][df['year'] == 1985]\n",
    "\n",
    "# Generate denominators of the weights\n",
    "\n",
    "df['g1_mean'] = df['g1'].groupby(df['year']).transform(lambda x: np.mean(x))\n",
    "\n",
    "df['g1_cont_1991mean'] = df['pg1_1991'].groupby(df['year']).transform(lambda x: np.mean(\n",
    "    (1 - df['g1'])*df['pg1_1991']/(1 - df['pg1_1991'])))\n",
    "\n",
    "# Get weights\n",
    "\n",
    "df['w1'] = df['g1']/df['g1_mean']\n",
    "\n",
    "df['w0'] = ((1 - df['g1'])*df['pg1_1991']/(1 - df['pg1_1991']))/df['g1_cont_1991mean']\n",
    "\n",
    "# Generate each component of the DID\n",
    "\n",
    "att_11 = np.mean(df['w1']*df['ypost'])\n",
    "att_10 = np.mean(df['w1']*df['ypre'])\n",
    "att_01 =  np.mean(df['w0']*df['ypost'])\n",
    "att_00 = np.mean(df['w0']*df['ypre'])\n",
    "\n",
    "# Get the ATT(1986,1986)\n",
    "att1986_1986 = (att_11 - att_10) -(att_01 - att_00)\n",
    "att1986_1986"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8953c0fe",
   "metadata": {},
   "source": [
    "* Calculate ATT(1986,1987)\n",
    "* This the the ATT for the group first treated at 1986 (g1==1) in the second period since treatment (1987)\n",
    "* Gen outcomes (easiest way to transform data into wide form is this)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 184,
   "id": "edad95ba",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "19.97935546103527"
      ]
     },
     "execution_count": 184,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df['ypost'] = df['y'][df['year'] == 1987] \n",
    "df['ypre'] = df['y'][df['year'] == 1985]\n",
    "\n",
    "df['g1_mean'] = df['g1'].groupby(df['year']).transform(lambda x: np.mean(x))\n",
    "df['g1_cont_1991mean'] = df['pg1_1991'].groupby(df['year']).transform(\n",
    "    lambda x: np.mean((1 - df['g1'])*df['pg1_1991']/(1 - df['pg1_1991'])))\n",
    "\n",
    "df['w1'] = df['g1']/df['g1_mean']\n",
    "df['w0'] = ((1 - df['g1'])*df['pg1_1991']/(1 - df['pg1_1991']))/df['g1_cont_1991mean']\n",
    "\n",
    "att_11 = np.mean(df['w1']*df['ypost'])\n",
    "att_10 = np.mean(df['w1']*df['ypre'])\n",
    "att_01 =  np.mean(df['w0']*df['ypost'])\n",
    "att_00 = np.mean(df['w0']*df['ypre'])\n",
    "\n",
    "att1986_1987 = (att_11 - att_10) -(att_01 - att_00)\n",
    "att1986_1987"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 185,
   "id": "90634f88",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "30.000872605079834"
      ]
     },
     "execution_count": 185,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df['ypost'] = df['y'][df['year'] == 1988] \n",
    "df['ypre'] = df['y'][df['year'] == 1985]\n",
    "\n",
    "df['g1_mean'] = df['g1'].groupby(df['year']).transform(lambda x: np.mean(x))\n",
    "df['g1_cont_1991mean'] = df['pg1_1991'].groupby(df['year']).transform(\n",
    "    lambda x: np.mean((1 - df['g1'])*df['pg1_1991']/(1 - df['pg1_1991'])))\n",
    "\n",
    "df['w1'] = df['g1']/df['g1_mean']\n",
    "df['w0'] = ((1 - df['g1'])*df['pg1_1991']/(1 - df['pg1_1991']))/df['g1_cont_1991mean']\n",
    "\n",
    "att_11 = np.mean(df['w1']*df['ypost'])\n",
    "att_10 = np.mean(df['w1']*df['ypre'])\n",
    "att_01 =  np.mean(df['w0']*df['ypost'])\n",
    "att_00 = np.mean(df['w0']*df['ypre'])\n",
    "\n",
    "att1986_1988 = (att_11 - att_10) -(att_01 - att_00)\n",
    "att1986_1988"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 186,
   "id": "ccf705c0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.636514\n",
      "         Iterations 4\n"
     ]
    }
   ],
   "source": [
    "logit_mod = sm.Logit.from_formula('g1 ~ 1', data = df.query(\n",
    "    '(year == 1997) & ((g1 == 1) | (time_til < 0))')).fit()\n",
    "df['pg1_1997'] = logit_mod.predict(df['g1']).g1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 187,
   "id": "aebe0bda",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "79.98798332042038"
      ]
     },
     "execution_count": 187,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df['ypost'] = df['y'][df['year'] == 1993] \n",
    "df['ypre'] = df['y'][df['year'] == 1985]\n",
    "\n",
    "df['g1_mean'] = df['g1'].groupby(df['year']).transform(lambda x: np.mean(x))\n",
    "df['g1_cont_1997mean'] = df['pg1_1997'].groupby(df['year']).transform(\n",
    "    lambda x: np.mean((1 - df['g1'])*(1 - df['g2'])*df['pg1_1997']/(1 - df['pg1_1997'])))\n",
    "\n",
    "df['w1'] = df['g1']/df['g1_mean']\n",
    "df['w0'] = ((1 - df['g1'])*(1 - df['g2'])*df['pg1_1997']/(1 - df['pg1_1997']))/df['g1_cont_1997mean']\n",
    "\n",
    "att_11 = np.mean(df['w1']*df['ypost'])\n",
    "att_10 = np.mean(df['w1']*df['ypre'])\n",
    "att_01 =  np.mean(df['w0']*df['ypost'])\n",
    "att_00 = np.mean(df['w0']*df['ypre'])\n",
    "\n",
    "att1986_1993 = (att_11 - att_10) -(att_01 - att_00)\n",
    "att1986_1993"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 188,
   "id": "19f21a73",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.693147\n",
      "         Iterations 1\n"
     ]
    }
   ],
   "source": [
    "logit_mod = sm.Logit.from_formula('g1 ~ 1', data = df.query('(year == 2003) & ((g1 == 1) | (time_til < 0))')).fit()\n",
    "df['pg1_2003'] = logit_mod.predict(df['g1']).g1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 189,
   "id": "ac684331",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "129.97868247120826"
      ]
     },
     "execution_count": 189,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df['ypost'] = df['y'][df['year'] == 1998] \n",
    "df['ypre'] = df['y'][df['year'] == 1985]\n",
    "\n",
    "df['g1_mean'] = df['g1'].groupby(df['year']).transform(lambda x: np.mean(x))\n",
    "df['g1_cont_2003mean'] = df['pg1_1997'].groupby(df['year']).transform(\n",
    "    lambda x: np.mean(df['g4']*df['pg1_2003']/(1 - df['pg1_2003'])))\n",
    "\n",
    "df['w1'] = df['g1']/df['g1_mean']\n",
    "df['w0'] = (df['g4']*df['pg1_2003']/(1 - df['pg1_2003']))/df['g1_cont_2003mean']\n",
    "\n",
    "att_11 = np.mean(df['w1']*df['ypost'])\n",
    "att_10 = np.mean(df['w1']*df['ypre'])\n",
    "att_01 =  np.mean(df['w0']*df['ypost'])\n",
    "att_00 = np.mean(df['w0']*df['ypre'])\n",
    "\n",
    "att1986_1998 = (att_11 - att_10) -(att_01 - att_00)\n",
    "att1986_1998"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 190,
   "id": "ab418588",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.636514\n",
      "         Iterations 4\n"
     ]
    }
   ],
   "source": [
    "logit_mod = sm.Logit.from_formula('g2 ~ 1', data = df.query('(year == 1997) & ((g2 == 1) | (time_til < 0))')).fit()\n",
    "df['pg2_1997'] = logit_mod.predict(df['g2']).g2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 191,
   "id": "3b3833d8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7.951605812652414"
      ]
     },
     "execution_count": 191,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df['ypost'] = df['y'][df['year'] == 1992] \n",
    "df['ypre'] = df['y'][df['year'] == 1991]\n",
    "\n",
    "df['g2_mean'] = df['g2'].groupby(df['year']).transform(lambda x: np.mean(x))\n",
    "df['g2_cont_1997mean'] = df['pg1_1997'].groupby(df['year']).transform(\n",
    "    lambda x: np.mean((1-df['g1'])*(1-df['g2'])*df['pg2_1997']/(1 - df['pg2_1997'])))\n",
    "\n",
    "df['w1'] = df['g2']/df['g2_mean']\n",
    "df['w0'] = ((1-df['g1'])*(1-df['g2'])*df['pg2_1997']/(1 - df['pg2_1997']))/df['g2_cont_1997mean']\n",
    "\n",
    "att_11 = np.mean(df['w1']*df['ypost'])\n",
    "att_10 = np.mean(df['w1']*df['ypre'])\n",
    "att_01 =  np.mean(df['w0']*df['ypost'])\n",
    "att_00 = np.mean(df['w0']*df['ypre'])\n",
    "\n",
    "att1992_1992 = (att_11 - att_10) -(att_01 - att_00)\n",
    "att1992_1992"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f5bcb49",
   "metadata": {},
   "source": [
    "### Jon Roth - Staggered"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 211,
   "id": "8a0853a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "\n",
    "utils.install_packages(\"staggered\")\n",
    "staggered = rpackages.importr('staggered')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 212,
   "id": "f1d9d18d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>estimate</th>\n",
       "      <th>se</th>\n",
       "      <th>se_neyman</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>68.328934</td>\n",
       "      <td>0.013269</td>\n",
       "      <td>0.013283</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    estimate        se  se_neyman\n",
       "1  68.328934  0.013269   0.013283"
      ]
     },
     "execution_count": 212,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "staggered_1 = staggered.staggered(df = rdf,\n",
    "                   i = 'id', g = 'treat_date', y = 'y', t = 'year', estimand = 'simple')\n",
    "\n",
    "with localconverter(ro.default_converter + pandas2ri.converter):\n",
    "    staggered_1 = ro.conversion.rpy2py(staggered_1)\n",
    "    \n",
    "staggered_1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 219,
   "id": "4e94c9a0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>estimate</th>\n",
       "      <th>se</th>\n",
       "      <th>se_neyman</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>8.657853</td>\n",
       "      <td>0.013073</td>\n",
       "      <td>0.013084</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   estimate        se  se_neyman\n",
       "1  8.657853  0.013073   0.013084"
      ]
     },
     "execution_count": 219,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "staggered_2 = staggered.staggered_cs(df = rdf,\n",
    "                   i = 'id', g = 'treat_date', y = 'y2', t = 'year', estimand = 'simple')\n",
    "\n",
    "with localconverter(ro.default_converter + pandas2ri.converter):\n",
    "    staggered_2 = ro.conversion.rpy2py(staggered_2)\n",
    "    \n",
    "staggered_2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 220,
   "id": "2a34e82d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>estimate</th>\n",
       "      <th>se</th>\n",
       "      <th>se_neyman</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>7.987738</td>\n",
       "      <td>0.013783</td>\n",
       "      <td>0.013804</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   estimate        se  se_neyman\n",
       "1  7.987738  0.013783   0.013804"
      ]
     },
     "execution_count": 220,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "staggered_3 = staggered.staggered_sa(df = rdf,\n",
    "                   i = 'id', g = 'treat_date', y = 'y2', t = 'year', estimand = 'cohort')\n",
    "\n",
    "with localconverter(ro.default_converter + pandas2ri.converter):\n",
    "    staggered_3 = ro.conversion.rpy2py(staggered_3)\n",
    "    \n",
    "staggered_3"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0fbb220b",
   "metadata": {},
   "source": [
    "* Aggregation over two ATTs(g,t)\n",
    "* Combine att1986_1986 att1986_1987\n",
    "* Formula is 2/T(T-1) Sum_g=2^T Sum_t=2^T 1(g<=t)ATT(g,t)\n",
    "*\tT is the number of g units (here 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 192,
   "id": "ce5fb152",
   "metadata": {},
   "outputs": [],
   "source": [
    "mean = 2/4*(4-1)\n",
    "sum1 = att1986_1986 + att1986_1987\n",
    "att = mean * sum1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da1c51b4",
   "metadata": {},
   "source": [
    "### Callaway & Sant'Anna via regressions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db621a45",
   "metadata": {},
   "source": [
    "* Calculate ATT(1986,1986)\n",
    "* This the the ATT for the group first treated at 1986 (g1==1) in the first period since treatment (1986)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 193,
   "id": "a29d63fc",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                          PanelOLS Estimation Summary                           \n",
      "================================================================================\n",
      "Dep. Variable:                      y   R-squared:                        0.9960\n",
      "Estimator:                   PanelOLS   R-squared (Between):              0.6936\n",
      "No. Observations:                2000   R-squared (Within):               0.9960\n",
      "Date:                Tue, Feb 15 2022   R-squared (Overall):              0.8541\n",
      "Time:                        22:11:01   Log-likelihood                    620.22\n",
      "Cov. Estimator:                Robust                                           \n",
      "                                        F-statistic:                    1.23e+05\n",
      "Entities:                        1000   P-value                           0.0000\n",
      "Avg Obs:                       2.0000   Distribution:                   F(2,998)\n",
      "Min Obs:                       2.0000                                           \n",
      "Max Obs:                       2.0000   F-statistic (robust):          1.343e+05\n",
      "                                        P-value                           0.0000\n",
      "Time periods:                       7   Distribution:                   F(2,998)\n",
      "Avg Obs:                       285.71                                           \n",
      "Min Obs:                       0.0000                                           \n",
      "Max Obs:                      1000.00                                           \n",
      "                                                                                \n",
      "                                Parameter Estimates                                \n",
      "===================================================================================\n",
      "                 Parameter  Std. Err.     T-stat    P-value    Lower CI    Upper CI\n",
      "-----------------------------------------------------------------------------------\n",
      "Intercept           8.4029     0.0079     1057.8     0.0000      8.3873      8.4185\n",
      "C(year)[T.1986]     1.0011     0.0132     76.081     0.0000      0.9753      1.0270\n",
      "treat               10.007     0.0252     397.38     0.0000      9.9574      10.056\n",
      "===================================================================================\n",
      "\n",
      "F-test for Poolability: 66.973\n",
      "P-value: 0.0000\n",
      "Distribution: F(999,998)\n",
      "\n",
      "Included effects: Entity\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\sajidmu2\\Anaconda3\\lib\\site-packages\\linearmodels\\panel\\model.py:1831: AbsorbingEffectWarning: \n",
      "Variables have been fully absorbed and have removed from the regression:\n",
      "\n",
      "g1\n",
      "\n",
      "  warnings.warn(\n"
     ]
    }
   ],
   "source": [
    "model = PanelOLS.from_formula(\"y ~ 1 + C(year) + g1 + treat + EntityEffects\", check_rank = False,\n",
    "                               drop_absorbed = True,\n",
    "                               data = df.query(\n",
    "                                   '((year==1985 | year==1986) & ((g1==1) | (time_til<0)))')).fit(cov_type = 'robust')\n",
    "print(model)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4753d5ed",
   "metadata": {},
   "source": [
    "* Now calculate ATT (1986, 1987)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 194,
   "id": "3716ba32",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                          PanelOLS Estimation Summary                           \n",
      "================================================================================\n",
      "Dep. Variable:                      y   R-squared:                        0.9991\n",
      "Estimator:                   PanelOLS   R-squared (Between):              0.8989\n",
      "No. Observations:                2000   R-squared (Within):               0.9991\n",
      "Date:                Tue, Feb 15 2022   R-squared (Overall):              0.9586\n",
      "Time:                        22:11:02   Log-likelihood                    699.77\n",
      "Cov. Estimator:                Robust                                           \n",
      "                                        F-statistic:                   5.316e+05\n",
      "Entities:                        1000   P-value                           0.0000\n",
      "Avg Obs:                       2.0000   Distribution:                   F(2,998)\n",
      "Min Obs:                       2.0000                                           \n",
      "Max Obs:                       2.0000   F-statistic (robust):          4.582e+05\n",
      "                                        P-value                           0.0000\n",
      "Time periods:                       8   Distribution:                   F(2,998)\n",
      "Avg Obs:                       250.00                                           \n",
      "Min Obs:                       0.0000                                           \n",
      "Max Obs:                      1000.00                                           \n",
      "                                                                                \n",
      "                                Parameter Estimates                                \n",
      "===================================================================================\n",
      "                 Parameter  Std. Err.     T-stat    P-value    Lower CI    Upper CI\n",
      "-----------------------------------------------------------------------------------\n",
      "Intercept           8.4029     0.0076     1100.7     0.0000      8.3880      8.4179\n",
      "C(year)[T.1987]     2.0107     0.0121     165.97     0.0000      1.9870      2.0345\n",
      "treat               19.979     0.0263     760.18     0.0000      19.928      20.031\n",
      "===================================================================================\n",
      "\n",
      "F-test for Poolability: 72.954\n",
      "P-value: 0.0000\n",
      "Distribution: F(999,998)\n",
      "\n",
      "Included effects: Entity\n"
     ]
    }
   ],
   "source": [
    "model = PanelOLS.from_formula(\"y ~ 1 + C(year) + g1 + treat + EntityEffects\", check_rank = False,\n",
    "                               drop_absorbed = True,\n",
    "                               data = df.query(\n",
    "                                   '((year==1985 | year==1987) & ((g1==1) | (time_til<0)))')).fit(cov_type = 'robust')\n",
    "print(model)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8b214b0",
   "metadata": {},
   "source": [
    "* Now calculate ATT (1986, 1988)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 195,
   "id": "c3ef04e1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                          PanelOLS Estimation Summary                           \n",
      "================================================================================\n",
      "Dep. Variable:                      y   R-squared:                        0.9995\n",
      "Estimator:                   PanelOLS   R-squared (Between):              0.9525\n",
      "No. Observations:                2000   R-squared (Within):               0.9995\n",
      "Date:                Tue, Feb 15 2022   R-squared (Overall):              0.9812\n",
      "Time:                        22:11:02   Log-likelihood                    614.16\n",
      "Cov. Estimator:                Robust                                           \n",
      "                                        F-statistic:                   1.099e+06\n",
      "Entities:                        1000   P-value                           0.0000\n",
      "Avg Obs:                       2.0000   Distribution:                   F(2,998)\n",
      "Min Obs:                       2.0000                                           \n",
      "Max Obs:                       2.0000   F-statistic (robust):          1.044e+06\n",
      "                                        P-value                           0.0000\n",
      "Time periods:                       9   Distribution:                   F(2,998)\n",
      "Avg Obs:                       222.22                                           \n",
      "Min Obs:                       0.0000                                           \n",
      "Max Obs:                      1000.00                                           \n",
      "                                                                                \n",
      "                                Parameter Estimates                                \n",
      "===================================================================================\n",
      "                 Parameter  Std. Err.     T-stat    P-value    Lower CI    Upper CI\n",
      "-----------------------------------------------------------------------------------\n",
      "Intercept           8.4029     0.0080     1054.6     0.0000      8.3873      8.4186\n",
      "C(year)[T.1988]     3.0015     0.0129     232.82     0.0000      2.9762      3.0268\n",
      "treat               30.001     0.0265     1132.3     0.0000      29.949      30.053\n",
      "===================================================================================\n",
      "\n",
      "F-test for Poolability: 66.744\n",
      "P-value: 0.0000\n",
      "Distribution: F(999,998)\n",
      "\n",
      "Included effects: Entity\n"
     ]
    }
   ],
   "source": [
    "model = PanelOLS.from_formula(\"y ~ 1 + C(year) + g1 + treat + EntityEffects\", check_rank = False, \n",
    "                              drop_absorbed = True,\n",
    "                               data = df.query(\n",
    "                                   '((year==1985 | year==1988) & ((g1==1) | (time_til<0)))')).fit(cov_type = 'robust')\n",
    "print(model)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ba86c3d",
   "metadata": {},
   "source": [
    "* Think now you got the idea!\n",
    "\n",
    "* Let me illustrate how this works for second group, ATT(1992,1992)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 196,
   "id": "27f52b22",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                          PanelOLS Estimation Summary                           \n",
      "================================================================================\n",
      "Dep. Variable:                      y   R-squared:                        0.9952\n",
      "Estimator:                   PanelOLS   R-squared (Between):              0.6199\n",
      "No. Observations:                1500   R-squared (Within):               0.9952\n",
      "Date:                Tue, Feb 15 2022   R-squared (Overall):              0.8263\n",
      "Time:                        22:11:03   Log-likelihood                    421.15\n",
      "Cov. Estimator:                Robust                                           \n",
      "                                        F-statistic:                   7.689e+04\n",
      "Entities:                        1000   P-value                           0.0000\n",
      "Avg Obs:                       1.5000   Distribution:                   F(2,748)\n",
      "Min Obs:                       0.0000                                           \n",
      "Max Obs:                       2.0000   F-statistic (robust):          8.521e+04\n",
      "                                        P-value                           0.0000\n",
      "Time periods:                      13   Distribution:                   F(2,748)\n",
      "Avg Obs:                       115.38                                           \n",
      "Min Obs:                       0.0000                                           \n",
      "Max Obs:                       750.00                                           \n",
      "                                                                                \n",
      "                                Parameter Estimates                                \n",
      "===================================================================================\n",
      "                 Parameter  Std. Err.     T-stat    P-value    Lower CI    Upper CI\n",
      "-----------------------------------------------------------------------------------\n",
      "Intercept           14.404     0.0094     1524.4     0.0000      14.386      14.423\n",
      "C(year)[T.1992]     1.0115     0.0168     60.303     0.0000      0.9786      1.0445\n",
      "treat               7.9516     0.0276     287.86     0.0000      7.8974      8.0058\n",
      "===================================================================================\n",
      "\n",
      "F-test for Poolability: 64.175\n",
      "P-value: 0.0000\n",
      "Distribution: F(749,748)\n",
      "\n",
      "Included effects: Entity\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\sajidmu2\\Anaconda3\\lib\\site-packages\\linearmodels\\panel\\model.py:1831: AbsorbingEffectWarning: \n",
      "Variables have been fully absorbed and have removed from the regression:\n",
      "\n",
      "g2\n",
      "\n",
      "  warnings.warn(\n"
     ]
    }
   ],
   "source": [
    "model = PanelOLS.from_formula(\"y ~ 1 + C(year) + g2 + treat + EntityEffects\", check_rank = False,\n",
    "                               drop_absorbed = True,\n",
    "                               data = df.query(\n",
    "                                   '((year== 1991 | year == 1992) & ((g2 == 1) | (time_til<0)))')).fit(cov_type = 'robust')\n",
    "print(model)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2900465f",
   "metadata": {},
   "source": [
    "* Now calculate  ATT(1992,1993)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 197,
   "id": "1230f97a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                          PanelOLS Estimation Summary                           \n",
      "================================================================================\n",
      "Dep. Variable:                      y   R-squared:                        0.9988\n",
      "Estimator:                   PanelOLS   R-squared (Between):              0.8695\n",
      "No. Observations:                1500   R-squared (Within):               0.9988\n",
      "Date:                Tue, Feb 15 2022   R-squared (Overall):              0.9509\n",
      "Time:                        22:11:03   Log-likelihood                    448.08\n",
      "Cov. Estimator:                Robust                                           \n",
      "                                        F-statistic:                   3.213e+05\n",
      "Entities:                        1000   P-value                           0.0000\n",
      "Avg Obs:                       1.5000   Distribution:                   F(2,748)\n",
      "Min Obs:                       0.0000                                           \n",
      "Max Obs:                       2.0000   F-statistic (robust):          3.288e+05\n",
      "                                        P-value                           0.0000\n",
      "Time periods:                      14   Distribution:                   F(2,748)\n",
      "Avg Obs:                       107.14                                           \n",
      "Min Obs:                       0.0000                                           \n",
      "Max Obs:                       750.00                                           \n",
      "                                                                                \n",
      "                                Parameter Estimates                                \n",
      "===================================================================================\n",
      "                 Parameter  Std. Err.     T-stat    P-value    Lower CI    Upper CI\n",
      "-----------------------------------------------------------------------------------\n",
      "Intercept           14.404     0.0093     1552.0     0.0000      14.386      14.422\n",
      "C(year)[T.1993]     2.0040     0.0162     123.93     0.0000      1.9722      2.0357\n",
      "treat               15.998     0.0277     577.99     0.0000      15.943      16.052\n",
      "===================================================================================\n",
      "\n",
      "F-test for Poolability: 66.027\n",
      "P-value: 0.0000\n",
      "Distribution: F(749,748)\n",
      "\n",
      "Included effects: Entity\n"
     ]
    }
   ],
   "source": [
    "model = PanelOLS.from_formula(\"y ~ 1 + C(year) + g2 + treat + EntityEffects\", check_rank = False,\n",
    "                               drop_absorbed = True,\n",
    "                               data = df.query(\n",
    "                                   '((year== 1991 | year == 1993) & ((g2 == 1) | (time_til<0)))')).fit(cov_type = 'robust')\n",
    "print(model)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb4af862",
   "metadata": {},
   "source": [
    "* For the third group, we ATT(1998,1998)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 198,
   "id": "94e4a17b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                          PanelOLS Estimation Summary                           \n",
      "================================================================================\n",
      "Dep. Variable:                      y   R-squared:                        0.9952\n",
      "Estimator:                   PanelOLS   R-squared (Between):              0.5053\n",
      "No. Observations:                1000   R-squared (Within):               0.9952\n",
      "Date:                Tue, Feb 15 2022   R-squared (Overall):              0.7932\n",
      "Time:                        22:11:03   Log-likelihood                    334.80\n",
      "Cov. Estimator:                Robust                                           \n",
      "                                        F-statistic:                   5.132e+04\n",
      "Entities:                        1000   P-value                           0.0000\n",
      "Avg Obs:                       1.0000   Distribution:                   F(2,498)\n",
      "Min Obs:                       0.0000                                           \n",
      "Max Obs:                       2.0000   F-statistic (robust):          4.812e+04\n",
      "                                        P-value                           0.0000\n",
      "Time periods:                      19   Distribution:                   F(2,498)\n",
      "Avg Obs:                       52.632                                           \n",
      "Min Obs:                       0.0000                                           \n",
      "Max Obs:                       500.00                                           \n",
      "                                                                                \n",
      "                                Parameter Estimates                                \n",
      "===================================================================================\n",
      "                 Parameter  Std. Err.     T-stat    P-value    Lower CI    Upper CI\n",
      "-----------------------------------------------------------------------------------\n",
      "Intercept           20.414     0.0110     1860.6     0.0000      20.392      20.436\n",
      "C(year)[T.1998]     1.0094     0.0212     47.694     0.0000      0.9678      1.0510\n",
      "treat               5.9477     0.0310     191.67     0.0000      5.8868      6.0087\n",
      "===================================================================================\n",
      "\n",
      "F-test for Poolability: 71.753\n",
      "P-value: 0.0000\n",
      "Distribution: F(499,498)\n",
      "\n",
      "Included effects: Entity\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\sajidmu2\\Anaconda3\\lib\\site-packages\\linearmodels\\panel\\model.py:1831: AbsorbingEffectWarning: \n",
      "Variables have been fully absorbed and have removed from the regression:\n",
      "\n",
      "g3\n",
      "\n",
      "  warnings.warn(\n"
     ]
    }
   ],
   "source": [
    "model = PanelOLS.from_formula(\"y ~ 1 + C(year) + g3 + treat + EntityEffects\", check_rank = False, \n",
    "                              drop_absorbed = True,\n",
    "                              data = df.query(\n",
    "                        '((year== 1997 | year == 1998) & ((g3 == 1) | (time_til<0)))')).fit(cov_type = 'robust')\n",
    "\n",
    "print(model)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d644647",
   "metadata": {},
   "source": [
    "* For the last group, we can't identify their ATT's because there is no valid comparison group!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
